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ABSTRACT 

We show that photoevaporation of small gaseous exoplanets (“mini-Neptunes”) in the habitable 
zones of M dwarfs can remove several Earth masses of hydrogen and helium from these planets and 
transform them into potentially habitable worlds. We couple X-ray/extreme ultraviolet (XUV)-driven 
escape, thermal evolution, tidal evolution and orbital migration to explore the types of systems that 
may harbor such “habitable evaporated cores” (HECs). We find that HECs are most likely to form 
from planets with ~ 1 M 0 solid cores with up to about 50% H/He by mass, though whether or not 
a given mini-Neptune forms a HEC is highly dependent on the early XUV evolution of the host star. 
As terrestrial planet formation around M dwarfs by accumulation of local material is likely to form 
planets that are small and dry, evaporation of small migrating mini-Neptunes could be one of the 
dominant formation mechanisms for volatile-rich Earths around these stars. 


1. INTRODUCTION 


Due to their small radii and low luminosities, M dwarfs 
currently offer the best opportunity for the detection of 
terrestrial planets in the habitable zone (HZ), the region 
around a star w here liquid water can exist on the sur¬ 


face o f a planet (Kasting et al. 1993 
20131. Upcoming missions such as t 


Kopparapu et al. 
le Transiting Exo¬ 
planet Survey Satellite (TESS) and the repurposed Ke¬ 
pler spacecraft (K2) will be capable of detecting po- 
tentiall y habitable Earths and super-Earth s around M 
dwarfs (Ricker et al.|2010 Howell et al.|2014). In particu¬ 
lar, the detection ot potentially habitable planets around 
M dwarfs is easier because the habita ble zones of these 
stars can be as close in as ~ 0.02 AU (Kopparapu et al. 


20131. However, such proximity implies that terrestrial 
planets forming within the HZs of low mass stars are 
likely to be small (< 0.3M m) and form dry (Raymond 


et al. 2007 Lissauer 20071. Moreover, M dwarfs are 


extrerhely active when young, bombarding their plan¬ 
ets with high energy radiation and bursts of relativis¬ 
tic particles during flaring events, which can erode their 
atmospheres and potentially s terilize the surface (Scalo 
et al.pOOT Segura et al.|2010 ). Strong tidal heating and 
orbital evolution may furth er impact the habita bility of 
planets around these stars (Barnes et al. 2013). Many 
planets formed in situ in the HZ ot iVI dwarfs may there¬ 
fore be uninhabitable. 

However, planets need not form and remain in place. It 
is now commonly accepted that both disk-driven migra¬ 
tion and planet-planet interactions can lead to substan¬ 
tial orbital changes, potentially bringing planets from 
outside the snow line (the region of the disk beyond which 
water and other volatiles condense into ices, facilitating 
the formation of massive planetary cores) to within the 
HZ dHayashi et al. 1985 Ida & Lin 2008a|b Ogihara 


&: lda||2l)Uy r Disk-driveii migration relies on the ex¬ 
change of angular momentum between a planet and the 


surrounding gaseous disk, which induces rapid i nward 
migra tion for planets in the terrestrial mass range (Ward 


1997). Around solar-type stars, dis k dispersal typically 


occurs on the orde r of a few Myr (Walter et al. 1988 
Strom et al. 1989). While there is evidence t hat disk 
lifetimes may exceed ^ 5 Myr for low mass stars (Carpen- 
20061 |Pascucci et al.||2009|) , disks are no longer 
present afte r ~ 10-20 Myr around all stellar types (Ribas 
et al.|[20l4 1. Planets migrating via interactions with the 


disk will thus settle into their new orbits relatively early. 
Planet-planet scattering, on the other hand, is by nature 
a stochastic process and may take place at any point dur¬ 
ing a system’s lifetime. However, such interactions are 
far more frequent during the final stages of planet assem¬ 
bly (up to a few tens of Myr), after whic h planets relax 
into their long-term quasi-stable orbits (Ford & Rasio 

2008t. 


Given the abundance of ices beyond the snow line, 
planets that migrate into the HZ from the outer regions 
of the disk should have abundant water, therefore satis¬ 
fying that important criterion for habitability within the 
HZ. However, the situation is complicated by the fact 
that these planets may have accumulated thick gaseous 
envelopes, which could render them uninhabitable. In¬ 
vestigating whether these so-called “mini-Neptunes” can 
lose their envelopes and form planets with solid surfaces 
is therefore critical to understanding the habitability of 
planets around low mass stars. 

In this study we focus on mini-Neptunes with initial 
masses in the range 1M 0 < Mp < 10 M 0 with up to 
50% hydrogen/helium by mass that have migrated into 
the HZs of mid- to late M dwarf stars. We investigate 
whether it is possible for atmospheric escape processes 
to remove the thick H/He envelopes of mini-Neptunes in 
the HZ, effectively turning them into volatile-rich Earths 
and super-Earths (terrestrial planets more massive than 
Earth),which we refer to as “habitable evaporated cores” 




















































































2 


Luger et al. 2015 


(HECs). We consider two atmospheric loss processes: 
XUV-driven escape, in which stellar X-ray/extreme ul¬ 
traviolet (XUV) photons heat the atmosphere and drive 
a hydrodynamic wind away from the planet, and a sim¬ 
ple model of Roche lobe overflow (RLO), in which the 
atmosphere is so extended that part of it lies exterior to 
the planet’s Roche lobe; that gas is therefore no longer 
gravitationally bound to the planet. We further couple 
the effects of atmospheric mass loss to the thermal and 
tidal evolution of these planets. Planets cool as they age, 
undergoing changes of up to an order of magnitude in ra¬ 
dius, which can greatly affect the mass loss rate. Tidal 
forces arising from the differential strength of gravity dis¬ 
tort both the planet and the star away from sphericity, 
introducing torques that lead to the evolution of the or¬ 
bital and spin parameters of both bodies. 


fects of atmospheric escape (Erkaev et al.||20071 |Murray- 

Clay et al.|20091|Tian|2009tlUwen & 

ackson|2U12| 

Lam- 

mer et al. 

2013D, Roche lobe overflow 

Trilling et al. 

1998 

Cu et al. 

2003), thermal evolution ( 

Lopez et al. 2UI2D, 

and tidal 

evolution (Jackson et al. 

SDDS'FiTTiTz- 

IVlello 

et al.||2008 Correia & Laskar||2011 ) c 

in exoplanets, 

none 


For some systems, in particular those that may harbor 
HECs, modeling the coupling of these processes is essen¬ 
tial to accurately determine the evolution, since several 
feedbacks can ensue. Tidal forces in the HZ typically 
act to decrease a planet’s semi-major axis, leading to 
higher stellar fluxes and faster mass loss. The mass loss, 
in turn, affects the rate of tidal evolution primarily via 
the changing planet radius, which is also governed by the 
cooling rate of the envelope. Changes to the star’s radius 
and luminosity lead to further couplings that need to be 
treated with care. 


Jackson et al. (2010) considered the effect of the cou- 
pling between rnass foss and tidal evolution on the hot 
super-Earth CoRoT-7 b, finding that the two effects are 
strongly linked and must be considered simultaneously 
to obtain an accurate understanding of the planet’s evo¬ 
lution. However, an analogous study has not been per¬ 
formed in the HZ, in great part because both tidal effects 
and atmospheric mass loss are generally orders of magni¬ 
tude weaker at such distances from the star. This is not 
necessarily true around M dwarfs, for two reasons: (a) 
their low luminosities result in a HZ that is much closer 
in, exposing planets to strong tidal effects and possible 
RLO; and (b) M dwarfs are extremely active early on, 
so that the XUV flux in the HZ can be orders of magni¬ 
tude hig her than that arou nd a solar-type star (see, for 
instance, Scalo et ^|2007 |. 

In this paper we present the results of the first model 
to couple tides, thermal evolution, and atmospheric mass 
loss in the habitable zone, showing that for certain sys¬ 
tems the coupling is key in determining the long-term 
evolution of the planet. We demonstrate that it is pos¬ 
sible to turn mini-Neptunes into HECs within the habit¬ 
able zone, providing an important pathway to the forma¬ 
tion of potentially habitable, volatile-rich planets around 
M dwarfs. 

In we provide a detailed description of the relevant 
physics. In we describe our model, followed by our 
results in ( ^ We then discuss our main findings and the 
corresponcling caveats in followed by a summary in 


f We present auxiliary derivations and calculations in 
3 Appendix. 

2. STELLAR AND PLANETARY EVOLUTION 

In the following sections we r eview the luminosity evo¬ 
lution of low mass stars ( ^2.1[ ), the habitable zone and 
its evolut ion (12.2), atmospheric escape processes from 
p lane ts (12.3), and tidal evolution of star-planet systems 
(«. 


2.1. Stellar Luminosity 

Following their formation from a giant molecular cloud, 
stars contract under their own gravity until they reach 
the main sequence, at which point the core temperature 
is high enough to ignite hydrogen fusion. While the Sun 
is thought to have spe nt < 50 Myr in th is pre-main 
sequence (PMS) phase (Baraffe et al. 1998), M dwarfs 
take hundreds of Myr to fully contract; the lowest mass 
M d warfs reach the main sequence only after ~ 1 Gyr 
(e.g., Reid fc Hawley|2005 ). During their contraction, M 
dwar is^emain at a rougEly constant effective tempera¬ 
ture (|Hayashi 1961), so that their luminosity is primar¬ 
ily a function of their surface area, which is significantly 
larger than when they arrive on the main sequence. M 
dwarfs therefore remain super-luminous for several hun¬ 
dred Myr, with total (bolometric) luminosities higher 
than the main sequence value by up to two orders of mag¬ 
nitude. As we discuss below, this will significantly affect 
the atmospheric evolution of any planets these stars may 
host. 

XUV emissions (lA A < lOOOA) from M dwarfs also 
vary strongly with time. This is because the XUV lu¬ 
minosity of M dwarfs is rooted in the vigorous ma gnetic 
fields gene rated in their large convection zones (Scalo 
et al. |20^ . Stellar magne tic fields are largely controlled 
by rotation (Parker 1955), whic h slows down wit h time 
due to angular momentum loss ( Skumanich||1972 ), lead¬ 
ing to an XUV activity that declines with stellar age. 
However, due to the difficulty of accurately determining 
both the XUV luminosities (usually inferred from line 
proxies) and the ages (often determined kinematically) 
of M dwarfs, the exact functional for m of the evolution 
is still very uncertain (for a review, see Scalo et al.|2007). 
Further complications arise due to the fact that the pro- 
cess(es) that generate magnetic fields in late M dwarfs 
may be quite different fr om those in e arlier type stars. 
The rotational dynamo of Parker] (1955) involves the am¬ 
plification of toroidal fields generated at the radiative- 
convective boundary; late M dwarfs, however, are fully 
convective, and have no such boundary. Instead, their 
magnetic fields may be formed by turbulent dynamos 
(Durney et al. 1993), which may evolve differently in 
time from those around higher mass stars (Reid & Haw- 
Ie^[200^ . 


Tn a comprehensive study of th e XUV emiss i ons of 
solar-type stars of different ages, Ribas et al. (2005) 
found that the XUV evolution is well modeled by a simple 
power law with an initial short-lived “saturation” phase: 


( -^/XUV j 

Lho\ 1 (_L_'^ 


(=)■ 


^ ^sat 
^ ^ ^sat5 


(1) 
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where Lxuv Lboi are the XUV and bolometric lumi¬ 
nosities, respectively, and (3 = —1.23. Prior to t = tsat, 
the XUV luminosity is said to be “saturated,” as obser¬ 
vations show that the ratio Uxuv !Lho\ remains relatively 
c onstant at early times. 

Jackson et al. (2012) find that 4at ~ 100 Myr and 
(J^xuv/^boi)sat ~ 10“^ for K dwarfs. Similar studies for 


M dwarfs, however, are still being developed (e.g., En¬ 
gle & Guinan 2011), but it is likely that the saturation 


times c ale is miuch longer for late-type M dwarfs. [Wright | 
et al. (2011|) show that X-ray emission in low mass stars 
IS saturated for Prot/T ^ 0.1, where Prot is the stellar ro¬ 
tation period and r is the convective turnover time. The 
extent of the convective zone increases with decreasing 
stella r mass; below 0.35Mfb, M dwarfs are fully convec¬ 
tive (Chabrier & B araffe 1997|), resulting in larger val¬ 
ues of r (see, e.g., Pizzolato 'et al. 200(1). Low mass 


stars also have longer spin-down tirnes (Stauffer et al. 
1994); together, these effects should lead to signihcantly 


longer saturation times compared to solar- type stars. 
This i s consistent with observational studies; [West et ah 
(2008) find that the magnetic activity lifetime increases 
from < 1 Gyr for early (i.e., most massive) M dwarfs 
to > 7 Gyr for late (least massive) M dwarfs, possibly 
due to the o nset of full convecti on around spectral type 
M3. Finally, Stelzer et al. (2013) find that for M dwarfs, 
/3 = —1.10 ±0.02 in the X-ray and /? = —0.79 ± 0.05 in 
the FUV (far ultraviolet), suggesting a slight ly shallower 
slope in the XUV compared to the value fromjRibas et al. 
([ 2 ^. 


2.2. The Habitable Zone 

The habitable zone (HZ) is classically defined as the 
region around a star where an Earth-like planet can sus¬ 
tain l iquid water on its surface (|Hart|1979 [Kasting et al. 
1993). Interior to the HZ, high surface teniperatures lead 
to the runaway evaporation of a planet’s oceans, which 
increases the atmospheric infrared opacity and reduces 
the ability of the surface to cool in a process known as 
the runaway greenhouse. Exterior to the HZ, greenhouse 
gases—gases, like water vapor, that absorb strongly in 
the infrared—are unable to maintain surface tempera¬ 
tures above the f reezing point, and the oceans freeze 
globally. Recently, Kopparapu et al. (2013) re-calculated 
the location of the HZ boundaries as a function of stel¬ 
lar luminosity and effective temperature using an up¬ 
dated one-dimensional, radiative-convective, cloud-free 
climate model. Their five boundaries are the (1) Recent 
Venus, (2) Runaway Greenhouse, (3) Moist Greenhouse, 
(4) Maximum Greenhouse, and (5) Early Mars habitable 
zones. 

The first and last limits can be considered “optimistic” 
empirical limits, since prior to ^ 1 and ^ 3.8 Gyr ago, re¬ 
spectively, Venus and Mars may have had liquid surface 
water. The ability of a planet to maintain liquid wa¬ 
ter and to sustain life at these limits is still unclear and 
probably depends on a host of properties of its climate. 
Conversely, the other three limits are the “pessimistic” 
theoretical limits, corresponding to where a cloud-free, 
Earth-like planet would lose its entire water inventory 
due to the greenhouse effect (2 and 3) and to where the 
addition of CO 2 to the atmosphere would be incapable 
of sustaining surface temperatures above freezing (4). 

Because stellar luminosities vary with time, the loca- 



a (AU) 


Fig. 1.— Location of the inner habitable zone (red), central hab¬ 
itable zone (green) and outer habitable zone (blue) as a function 
of stellar mass at fO Myr (dashed) and 1 Gyr (solid). After 1 Gyr, 
the evolution of the HZ is negligible for M dwarfs. 


tion of the HZ is not fixed. In Figure we plot the HZ 
at 10 Myr (dashed lines) and a t 1 Gyr (solid line s ), cal- 
culated from the HZ model of Ko pparapu et al. ([2013 ) 
and the stellar evolution models of 


Haratte et al. 


9981 


While for K and G dwarfs the change in the HZ is negligi¬ 
ble, the HZ of M dwarfs moves in substantially, changing 
by nearly an order of magnitude for the least massive 
stars. Due to this evolution, planets observed in the HZ 
of M dwarfs today likely spent a significant amount of 
time interior to the inner edge of the HZ, provided they 
either formed in situ or migrated to their cu rrent posi¬ 
tions relatively early, huger & Barnes (|2015 ) explore in 
detail the effects of the evolution of the HZ on terrestrial 
planets. 

Finally, we note that the location of the HZ is also a 
function of the eccentricity e. This is due to the fact that 
at a fixed semi-major axis a, the or bit-averaged flux (F) 


is hig her for more eccentric orbits (Williams & Pollard 


2002 ): 


(F) = 


Tbol 


4710^1/1 ~ 


( 2 ) 


2.3. Atmospheric Mass Loss 


Planetary atmospheres constantly evolve. Several 
mechanisms exist through which planets can lose signifi¬ 
cant quantities of their atmospheres to space, leading to 
dramatic changes in composition and in some cases com¬ 
plete atmospheric erosion. The early Earth itself could 
have been rich in hydrogen, with mix ing ratios as high 


as 30% in the prebioti c atmosphere (jTian et al. 2005 
Hashimoto et al. 2007). A variety of processes subse- 


gueni 


et al. 


iv led to the loss of most of this hydrogen; jWatsonj 


(1981) argue that on the order of 10^^ g of hydrogen 


could have escaped in the first billion years. Similarly, 


Kasting & Pollack (1983) calculated that early Venus 
could have lost an Earth ocean equivalent of water in the 
same amount of time. Currently, observational evidence 
for atmos pheric escape exists for two “hot Jupiters,” HD 
209458b (jVidal-Madjar et al.j [20031) and HD 189733b 
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(Lecavelier Des Etangs et al.||2010 ), and one “hot Nep- 
tune,” GJ 43tjb (Kulow et al.|2U14l, whose proximity to 
their parent stars leads to the rapid hydrodynamic loss 
of hydrogen. 

Atmospheric escape mechanisms fall into two major 
categories: thermal escape, in which the heating of 
the upper atmosphere accelerates the gas to velocities 
exceeding the escape velocity, and nonthermal escape, 
which encompasses a variety of mechanisms and may 
involve energy exchange among ions or erosion due to 
stellar winds. While nonthermal processes certainly do 
play a role in the evaporation of super-Earth and mini- 
Neptune atmospheres, the high exospheric temperatures 
resulting from strong XUV irradiation probably make 
thermal escape the dominant mechanism for planets 
around M dwarfs at early times. However, the escape 
can be greatly enhanced by flares and coronal mass ejec¬ 
tions, which can completely erode the a tmosphere of a 
planet lacking a stron g magnetic field (Lammer et al. 
20^ [Scab etar]|2007l ). For a review of the nonthermal 


mechanisms of escape, the reader is referred to Hunten 


2 . 3 . 1 . Jeans Escape 

In the low temperature limit, atmospheric mass loss 
occurs via the ballistic escape of individual atoms from 
the collisionless exosphere, where the low densities ensure 
that atoms with outward velocities exceeding the escape 
velocity will escape the planet. Originally developed by 
Jeans] (1925), who assumed a hydrostatic, isothermal at- 
he 


mosphere, the mas s loss rate of a pure hydrogen atmo¬ 
sphere is given by ( 6pik|[l963 1 


dMp 

dt 


= ATTRz^^nmHVt 


(1 + ^j)e 

2y^ 


— Aj 


( 3 ) 


where J?cxo is the radius of the exobase, n is the number 
density of hydrogen atoms at the exobase, mn is the 
mass of a hydrogen atom, vt is the thermal velocity of 
the gas, and Aj is the Jeans escape parameter, defined 
as the ratio of the potential energy to the kinetic energy 
of the gas and given by 


A,/ = 


GMpmH 


hT R 

exo-^ ^Gxo 


( 4 ) 


where G is the gravitational constant, Mp is the mass of 
the planet, and T^xo is the temperature in the (isother¬ 
mal) exosphere. Since in the Jeans regime the thermal 
velocity of the gas is less than the escape velocity, the 
bulk of the gas remains bound to the planet, and only 
atoms in the tail of the Maxwell-Boltzmann distribution 
are able to escape. Jeans escape is thus slow. As an 
example, the present Jeans escape flux f or hydrogen on 


Venus is on the order of 10 cm ^s ^ (jLammer et al. 


20061, corresponding to the feeble rate of ^ 10 g/s, 


which would remove only one part in 10^^ of the atmo¬ 
sphere per billion years. 

For higher exospheric temperatures and/or larger val¬ 
ues of i?exo! however, corresponding to low values of Aj, 
the atmosphere enters a regime in which the velocity of 
the average atom in the exosphere approaches the escape 
velocity of the planet. In this regime, the bulk of the 
upper atmosphere ceases to be hydrostatic (and isother¬ 


mal), Q is no longer valid, and the escape rate must be 
calculated from hydrodynamic models. 

2 . 3 . 2 . Hydrodynamic Escape 

One of the primary mechanisms for heating the exo¬ 
sphere and decreasing Aj is via strong XUV irradiation. 
XUV photons are absorbed close to the base of the ther¬ 
mosphere, where they deposit energy and heat the gas 
via the ionization of atomic hydrogen. This heating is 
balanced by the adiabatic expansion of the upper atmo¬ 
sphere, downward heat conduction, and cooling by re¬ 
combination radiation. For sufficiently high XUV fluxes, 
the expansion of the atmosphere accelerates the gas to 
supersonic speeds, at which point a hydrodynam ic wind 
is est ablished similar to the solar Parker wind (Parker] 
19651. Once the gas reaches the exosphere, it will escape 
the planet provided its kinetic energy exceeds the energy 
required to lift it out of the planet’s gravitational well. 

Since the kinetic energy of a hydrogen atom at the 
exobase is |fcTexo, the condition Xj < 1.5 implies that 
the kinetic energy of the gas is greater than the absolute 
value of its gravitational binding energy, and it should 
therefore begin to escape in bulk in a process commonly 
referred to as “blow-off.” Unlike in the Jeans escape 
regime, where the mass loss occurs on a per-particle ba¬ 
sis, blow-off leads to the rapid loss of large portions of 
the upper atmosphere, irrespective of particle species, as 
atoms and molecules heavier than hydrogen are carried 
along by the hydr odynamic wind. However, contrary to 
what Opik (19631 suggests, the mass loss in this stage is 


not arbitrarily high, since once blow-off begins the up¬ 
per atmosphere can no longer be treated as isothermal. 
As the exosphere expands it also cools, so that in the 
absence of an energy source the value of Aj will tend to 
increase, thereby moderating the blow-off. The mass loss 
is, in this sense, “energy-limited,” and may be calculated 
by equating the energy input to the energy required to 
drive the escape. 

Originally proposed by j Watson et al.j (|1981|), the 
energy-limited model assunies that the XUV hux is ab¬ 
sorbed in a thin layer at radius i?xuv where the optical 
depth to stellar XUV ph otons is unity. Recen tly updated 
to include tidal effects by Erkaev et al. (20071, this model 


approximates the mass loss as 


dMp ^ cxuvttFxuv^p^xuv 
dt 


GMpiftide(C) 


( 5 ) 


where exuv is the heating efficiency parameter (see be¬ 
low) , Exuv is the incident XUV flux, Rp is the planetary 
radius, and Eftide is a tidal enhancement factor, account¬ 
ing for the fact that for sufficiently close-in planets, the 
stellar gravity reduces the gravitational binding energy 
of the gas such that i t need only reach th e Roche radius 
to escape the planet. Erkaev et al. (20071 show that 


iftide(C) = 1 - + 


1 

2^ ' 2^ 

where the parameter ^ is defined as 

J^Roche 




Rxuv 


( 6 ) 


(7) 
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with 


R = ' 

-^Roche — 


3M* 


1/3 


a, 


( 8 ) 


where M* is the mass of the star and a is t he sem i-major 
axis. For simplicity, as in Lopez et al. (2012), we re¬ 


place Rp with i?xuv in ([^ j which is approximately valid 


Limu .tt.xuv 
Murray-Glay et al. 

2009 

Lopez et al. 

2012) 

Since the input 2 

GJV energy is ba 

anced 


cooling radiation (via Lyman a emission in the case 
of hydrogen) and heat conduction, only a fraction of 
it goes into the adiabatic expansion that drives escape. 
Rather than running complex hydrodynamic and radia- 


Lopezetal 

^ choose 


20 


se to folc 


Leitzinger et al. 

201l| 

2012 Lammer 

et al. 


and cooling into an efficiency parameter, exuv, defined 
as the fraction of the incoming XUV energy that is con¬ 
verted into PdV work. Because of the complex depen¬ 
dence of exuv on the atmospheric compositi on and struc¬ 
ture, its value is still poorly constrained. Lopez et al. 
(2012) estimate exuv = 0.2 ± 0.1 for super-Earths and 


mini-Neptunes based on va lues found through out the lit¬ 
erature. Earlier work by Chassefiere (1996) estimates 
exuv ^ 0.30 for Venus-like planets but tne author argues 
that the actual v a lue m ay be closer to 0.15. Recently, 
Owen & Jackson (2012) found X-ray efficiencies < 0.1 
tor planets more massive than Neptune, but highe r effi- 
ciencies (~ 0.15) for terrestrial planets. Moreover, She- 


matovich et al. (2014) argue that studies that assume et- 
hciencies higher than about 0.2 probably lead to overesti¬ 
mates in the escape rate. On the other hand, some stud¬ 
ies su ggest higher heating efficiencies: Koskinen et al. 
(2012) use hydrodynamic and photochemical models of 
the hot Jupiter HD209458b to calculate exuv = 0.44. 

As we have already implied , unlike Jeans esca pe, hy¬ 
drodynamic blow-off is fast. Chassefiere (1996) calcu¬ 
lates the maximum hydrodynamic escape rate from the 
early Venusian atmosphere to be ~ 10® g/s, ten orders 
of ma gnitud e higher than the present Jeans escape flux 
(see § 2.3.1). Although there has been debate over the 
validity ot th e blow-off assumpt ion (see, f or instance, the 
discus sion in Xian et al. 2008), recently Lammer et al. 
(2013) showed that for super-Earths exposed to high lev¬ 
els ot XUV irradiation, the energy-limited approximation 
yields mass loss rates that are consistent with hydrody¬ 
namic models to within a factor of about two, which is 
within the uncertainties of the problem. 

2 . 3 . 3 . Controlled Hydrodynamic Escape 

It is also worth noting that there may be an interme¬ 
diate regime between Jeans escape and blow-off known 
as “modi fied Jeans escape” o r “controlled hydrodynamic 
escape” (Erkaev et al. 2013). In this regime, which oc¬ 
curs for intermediate XUV fluxes and/or higher plan¬ 
etary surface gravity, blow-off conditions are not met 
but the atmosphere still expands, so that the hydrostatic 
Jeans formalism is not valid. In order to calculate the 
escape rate, one must replace the classical Maxwellian 
velocity distribution with one that includes the bulk ex¬ 
pansion velocity of the atmosphere. This yields escape 


rates lower than those due to a hydrodynamic flow, but 
significantly higher than those predicted by the hydro¬ 
static Jeans equation 

2 . 3 . 4 . Jeans Escape or Hydrodynamic Escape? 

Since the location of the habitable zone is governed 
primarily by the total (bolometric) flux incident on a 
planet, the higher ratio of Lxuv to Lboi of M dwarfs 
implies a much larger XUV flux in the HZ compared to 
solar-type stars. The present-day solar XU V luminosity 
is Lx uv/Lboi ~ 3.4 x 10“® (see Table 4 in Ribas et al.f 
2005), while for ac tive M dwarfs this ratio is ~ 10“'^ (e.g., 
Scalo et al. 2007). Therefore we should expect planets 
in the HZ of M dwarfs to experience XUV fluxes sev¬ 
eral orders of magnitude greater than the present Earth 
level (Exu vm ^ 4.64erg/s/cm^). Recent papers ( Lam-| 


mer et al. 2007[ 2013[ |Erkaev et al. 


2013) show that 


terrestrial planets experiencing XUV fluxes correspond¬ 
ing to 10 X and 100 x Exuv© are in the hydrodynamic 
flow regime, and we may thus expect the same for super- 
Earths/mini-Neptunes in the HZ of active M dwarfs. 

In Figure we plot the evolution of the XUV flux re¬ 
ceived by a panet located close to the inner edge of the 
HZ (defined at 5 Gyr), for three different M dw arf m asses 
and two different XUV saturation times (see (12.1|). The 


dashe d lines correspond to the critical fluxes m Erkaev 


et al. (2013) above which hydrodynamic escape occurs, 
tor 1 and 10 M 0 and two values of exuv • Earth-mass 
planets remain in the hydrodynamic escape regime for 
at least 1 Gyr in all cases and in excess of 10 Gyr for 
active M dwarfs. The duration of this regime is shorter 
for 10 Mq planets, but still on the order of several 100 
Myr. 

We note also that tidal effects can significantly increase 
the critical value of Aj be low which hydrodyn amic escape 
occurs (see discussion in Erkaev et al. 2007). Hydrody¬ 
namic escape ensues wheh the thermaf energy of the gas 
exceeds its potential energy, which occurs when 


1.5 > 


GAJpKiidQTTipf 

ETexo77exo 


or 


AJ < 


1.5 


Eftide 
= Acrit 


(9) 


provided we maintain the original definition of the Jeans 
parameter Due to the strong tidal forces acting on 
the planets we consider here, this effect should greatly 
increase the critical value of the escape parameter, ef¬ 
fectively reducing the value of Exuv for which hydrody¬ 
namic escape occurs. 

2 . 3 . 5 . Energy-Limited or 
Radiation/Recombination-Limited? 

Hydrodynamic escape from planetary atmospheres 
need not be energy-limited. In the limit of high extreme 
ultraviolet (EUV) flu x (low-energy XUV photo ns with 
lOOA < A < lOOOA), |Murray-Clay et al. (2009) showed 
that the escape is “radiation/recombination-limited,” 

scaling roughly as M oc (Eeuv)^^^- In this regime, the 
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Fig. 2.— Evolution of the XUV flux received by planets close to the inner edge of the HZ (at 1 Gyr) for stars of mass 0.1, 0.2, and 0.3 
Mq. Solid lines correspond to an XUV saturation time of 0.1 Gyr; dashed lines correspond to 1 Gyr. The flux at Earth is indicated by the 
black line. The dot corresponds to the earliest time for which [Ribas et ar] l |200^ has data for solar-type stars; this is also roughly the time 
at which Earth formed. An XUV luminosity saturated at 10 '^Lboi is roughly indicated by the dotted line. Finally, the dashed gray lines 
indicate the minimum XUV fluxes required to sustain blowoff according to the study ofjErkaev et al.j (|2013||. Super-Earths in the HZs of 
M dwarfs remain in the blowoff regime for at least a few 100 Myr; Earths undergo blowott tor mucii longer. For an XUV saturation time 
of 1 Gyr, blowoff occurs for several Gyr for all planets. 


upper atmosphere thermostats to T ^ 10"^ K, photoion¬ 
ization is balanced by radiative recombination (as op¬ 
posed to PdV work), and a large fraction of the gas en¬ 
ergy budget is lost to Lyman a cooling. The mass loss 
rate is found by solving the mass continuity equation, 
yielding 

M = inr^psCs (10) 


where Vs is the radius of the sonic point (where the 
wind velocity equals the local sound speed c^) and Ps 
is the density at Vg- Since the bulk of the flow is ion¬ 
ized, the density is fixed by ionization-recombination bal- 

ance, scaling roughly as (T’euv^'s) ■ For a 0.7Mj hot 
Jupiter , the radiation/recombin ation limited mass loss 
rate is ( Murray-Clay et ar]|2009 ) 


M, 


RR 


4 X 10^^ g/s 


F, 


EUV 


5 X 10® erg/cm^/s 


1/2 


( 11 ) 


Owen & Jackson (2012) re-derive this expression with 
explicit scalings on several planet properties: 


Mrr « 2.4 X 10“ g/s (1 -b /parker 






Rr, 


10 Rf, 


3/2 / ^ \ 1/2 

1/3/ ly 10 km/s 


1040 s-i 
ceuv 


1/2 


(12) 


where the quantity (1 -I- x) is the radius of the ioniza¬ 
tion front in units of i?p, /parker is the Mach number of 
the flow, is the stellar EUV luminosity in photons per 
second, A is a geometrical factor and ceuv is the isother¬ 
mal sound speed of the gas. Taking x = 0, /parker = 1, 
A = 1/3, Ceuv = 10 km/s, and an average EUV photon 
energy hv = 20 eV, this becomes 


Mrr « 7.11 X 10^ g/s 


Feuv 

erg/cm^/s 


1/2 


Rff 


3/2 


(13) 


The transition from energy-limited to radiation/ 
recombination-limited esc^e is fo und by equating the 
two expressions, namely and (13), and solving for 
the critical EUV flux. Tor hot Jupiters and mini- 
Neptunes alike, the transition occurs at roughly Feuv 
104 erg/cm^/s. Below this flux, the escape is energy- 
limited; above it, the escape is radiation/recombination- 
limited and thus increases more slowly with the flux. 

Mini-Neptunes that migrate early into the HZ of M 
dwarfs are exposed to EUV fluxes up to an order of mag¬ 
nitude larger than this critical value. During this period, 
which lasts on the order of a few hundred Myr, the mass 
loss rate may be radiation/recombination-limited. 

However, whether the hi^-flux mas s loss rate is more 
accurately described by ® or (13) will depend on 
wheth er the flux is dominated by iGray or EUV radi¬ 
ation. Owen & Jackson (2012) show that the mass loss 
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rate for X-ray driven hydrodynamic winds scales linearly 
with the X-ray flux; this is because the sonic point for 
X-ray flows tends to occur below the ionization front. 
Even though recombination radiation still removes en¬ 
ergy from the flow, it does so once the gas is already 
supersonic and thus causally decoupled from the planet, 
such that it cannot bottleneck the escape. Although the 
authors caution that the dependence of the mass loss 
rate on planet mass and radius must be determined nu¬ 
merically, this regime is analogous to the ener^-limited 
regime and can be roughly approximated by Q. 

For X-ray luminosities Lx ^ IO^^Lq, close-in plan¬ 
ets under go X-ray driven hydrody namic escape (see Fig¬ 
ure 11 in Owen & Jackson 2012). If X-rays contribute 
significantly to the XUV emissions of young M dwarfs, 
their X-ray luminosities may exceed this value for as 
long as 1 Gyr, and close-in mini-Neptunes will undergo 
energy-limited escape. Unfortunately, given the lack of 
observational constraints on the X-ray/EUV luminosi¬ 
ties of young M dwarfs, it is unclear at this point whether 
the hydrodynamic escape will be EUV-driven (radiation/ 
recombination-limited) or X-ray-driven (energy-limited). 

2 . 3 . 6 . The Effect of Eecentricity 

Most of the formalism that has been developed to ana¬ 
lytically treat hydrodynamic blow-off only considers cir¬ 
cular orbits. For planets on sufficiently eccentric orbits, 
neither the stellar flux nor the tidal effects may be treated 
as constant over the course of an orbit. Due to this fact, 
an expression analogous to ([^ for eccentric orbits seems 
to be lacking in the literature. In this section we derive 
such an expression. 

There are two separate effects that enhance the mass 
loss for planets on eccentric orbits. Most papers account 
for the first effect, which is the increase of the orbit- 
averaged stella£^iD^_by_a facto of l/Vl — (see, for 
instance, Kopparapu et al.||2013 ). However, for e < 0.3, 
this effect is quite small. T'he second, more important 
effect is that the Roche lobe radius is no longer constant 
over the course of an orbit, and ^ is not valid. Instead, 
we must replace a with the instantaneous planet-star sep¬ 
aration r(t): 


RRoche(0 — 


/ Mp 


V3M* 


1/3 


r{t). 


(14) 


One might wonder whether this replacement is valid. 
Specifically, if i?Roche(0 changes faster than the atmo¬ 
sphere is able to respond to the changes in the gravita¬ 
tional potential, we would expect that the time depen¬ 
dence of the mass loss rate would be a complicated func¬ 
tion of the tides generated in the atmosphere. On the 
other hand, if the orbital period is very large compared 
to the dynamical timescale of the planet, the atmosphere 
will have sufficient time to assume the equilibrium shape 
dictated by the new potential. This limit is know n as 
the quasi-static approximation ( Sepinsky et al.|2007 ). In 
the Appendix, we show that all the planets in our runs 
with ecce ntric ities e 0.4 are in the quasi-static regime 
and that (14) is therefore valid. In some runs, we allow 
the eccentricity to incr ease beyond 0.4. We discuss the 
implications of this in §5.9| 

Since RRoche — RRoche {t ), A^tidej and dilL/dt (Equa¬ 

tions 1^ and are now also functions of t, varying 


significantly over a single orbit. To account for this, we 
may calculate the time-averaged mass loss rate over the 
course of one orbit, {M)t, such that the total amount of 
mass lost in time At is AM = {M)tAt. To this end, in 
the Appendix we derive the eccentric version of ATtide: 


1 






r2TT 


(1 ecosE) 2^ ^ 253 ( 1 - 6 cosL;)2 

(15) 

where £ is the time-independent parameter given by ([^ 
and (1^ and E is the eccentric anomaly. The average 
mass loss rate is then simply 


(M)t = 


Mq 

K 

j- P CX 


where 


Mq = 


_ -^xuv^xuyAxuv 


AGMpO? 


(16) 


(17) 


is the zero-eccentricity mass loss rate in the absence of 
tidal en hancem ent. Note that the flux-enhancement fac¬ 
tor l/Vl — e2 is already folded into_Aecc, since it must 
be incorporated when integrating (A7). 


For 5 > 10, the integral may be approximated by the 
analytic expression 





Co 452 


— 


(18) 


which greatly reduces computing time. 

As we show in the Appendix, the decreased Roche 
lobe distance for eccentric orbits has a large effect on 
the amount of mass lost, pa rtic ularly for low values of 
5 and for high e (see Figure 13). Moreover, higher ec¬ 


centricities result in Roche lobe overflow at larger values 
of a compared to the circular case, since the planet may 
overflow near pericenter, leading to mass loss rates po¬ 
tentially orders of magnitude higher. We disc uss RLO in 
detail in the Model Description section ((3.3). 


2 . 4 . Tidal Theory 

The final aspect of planetary evolution we discuss is the 
effect of tidal interactions with the host star. Below we 
review two different approaches to analytically calculate 
the orbital evolution of the system. 

2 . 4 . 1 . Constant Phase Lag 

Classical tidal theory predicts that torques arising from 
interactions between tidal deformations on a planet and 
its host star lead to the secular evolution of the orbit and 
the spin of both bodies. In this paper we employ the 


“equilibrium tide” model of Darwin (1880), which ap¬ 
proximates the tidal potential as a superposition of Leg¬ 
endre polynomials representing waves propagating along 
the surfaces of the bodies; these add up to give the fa¬ 
miliar tidal “bulge.” Because of viscous forces in the 
bodies’ interiors, the tidal bulges do not instantaneously 
align with the line connecting the two bodies. Instead, 
the wave on the body will lag or lead by an angle 
eAT^i, assumed to be constant in the constant phase lag 
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(CPL) model. In general, different waves may have dif¬ 
ferent lag angles, and it is unclear how the eN,i v ary as a 


function of freque ncy. A common approach (see Ferraz- 


20081 is to assu me that the magnitudes o t 


es are equal (see Goldreich k. Soter 1966), 


Mello et al. 
the lag ang. 

while their signs may change depending on the orbital 
and rotational frequencies involved. This allows us to 
introduce the tidal quality factor 


Qi = —, 


(19) 


which in turn allows us to express the lags (in radians) 
as 

eN,^ = ±^. ( 20 ) 

The parameter Qi is a measure of the dissipation within 
the body; it is inversely proportional to the amount 
of orbital and rotational energy lost to heat per cycle, in 
analogy with a damped-driven harmonic oscillator. The 
merit of this approach is that the tidal response of a body 
can be captured in a single parameter. Planets with high 
values of Qp have smaller phase lags, dissipate less en¬ 
ergy and undergo slower orbital evolution; planets with 
low values of Qp have larger phase lags, higher dissipa¬ 
tion rates, and therefore faster evolution. Measurements 
in the solar system constrain the value of Qp for ter¬ 
restrial bodies in the range 10- 500, while gas giants ar e 
consistent with Qp ~ 10^ — 10® (Goldreich & Soter|l966 ). 
Values of Q* for the Sun and other main sequence stars 
are poorly constrained but are likely to be > 10® — 10® 
(Schlaufman et al. 2010| Penev et al.||2012 1. Intuitively, 
this makes sense, given that the dissipation due to inter¬ 
nal friction in rocky bodies should be much higher than 
that in bodies dominated by gaseous atmospheres. One 
should bear in mind, however, that the exact dependence 
of Qi on the properties of a body’s interior is likely to 
be extremely complicated. Given the dearth of data on 
the composition and internal structure of exoplanets, it 
is at this point impossible to infer precise values of Qp 
for these planets. 

By calculating the forces and torques due to the tides 
raised on both the planet and the star, one can arrive at 
the secular expressions for the evolution of the planet’s 
orbital parameters, which are given by a set of coupled 
nonlinear differential equations; these are reproduced in 
the Appendix. 

2.4.2. Constant Time Lag 

Unlike the GPL model, which assumes the phase lag of 
the tidal bulge is constant, the constant time lag (CTL) 
model assumes that it is the time interval between the 
bulge and the passage of the pe rturbing body th at is 
constant. Originally proposed b y |Alexander| (fl973| and 
updated by Leconte et al. (2010), this model allows for a 
continuum ot tidal wave frequencies and therefore avoids 
unphysical discontinuities present in the GPL model. 
However, implicit in the GTL theory is the assumption 
that the lag an gles are directly proportional to the driv¬ 
ing frequency (Greenberg 2009), which is likely also an 
oversimplification. We note, however, that in the low 
eccentricity limit, both the CPL and the CTL models 
arrive at qualitatively similar results. At higher eccen¬ 
tricities, the CTL model is probably better suited, given 


TABLE 1 

Free Parameters and Their Ranges 


Parameter 

Range 

Default 

Notes 



0.08 - 0.4 

_ 

Late-mid MD 

Mp{ Me) 

1-10 


- 



RxvviRp) 

1.0 - 1.2 

1.2 

See q2.3.2| 

a 

IHZ - OHZ 


See^XT 


e 

0.0 - 0.95 


- 



Po,* (days) 

1.0 - 100 

30.0 

Initial rot. per. 

Ih 

10-® - 0.5 


H mass fraction 

exuv 

0.1 - 0.4 

0.3 

- 



^min 

1 -1- 10-5 - 3 

3 

See [ 

3.3 


Atmos, esc. 

R/R-Lim / E-Lim 

- 

See f 

3.3 


Tidal model 

CPL/CTL 

CTL 

- 



Q* 

10® - 10® 

10® 

CPL only 

Qp 

lO^ - 10® 

lO^ 

CPL only 

r* (s) 

10-2 - 10-1 

10-1 

CTL only 

Tp (s) 

10-® - 10® 

10-1 

CTL onR 


/3 

0.7- 1.23 

1.23 

See Eq. dl 

] 

Isat (Gyr) 

0.1 - 1.0 

1.0 

XUV sat. mne 

to (Myr) 

10.0 - 100.0 

10.0 

Integration start 

istop (Gyr) 

0.01 - 5.0 

5.0 

Integration end 


that it is derived to eighth order in e (versus second order 
in the CPL model). 

The tidal quality factors Qi do not enter the CTL cal¬ 
culations at any point; instead, the dissipation is charac¬ 
terized by the time lags r^. Al though there is no g eneral 


conversion between Qi and r^, Leconte et al. (2010) show 


ution. 


( 21 ) 


that provided annual tides dominate the evo. 

1 

where n is the mean motion (or the orbital frequency) of 
the secondary body (in this case, the planet). 

For a planet with Qp = 10^ in the center of the HZ 
of a late M dwarf, Tp « 10 s; rocky planets with lower 
Qp may have values on the order of hundreds of seconds. 
-1 


time lags. For reference, [Leconte et al.| ( |2010[ ) argue that 
hot Jupiters should have 2 x 10“® s < Tp < 2 x 10“^ s. 

The tidal evolution expressions are reproduced in the 
Appendix. For a more detailed review of tidal theory. 


the reader is referred to Ferraz-Mello et al. _([2008|, 


et al. (2011), and the Appendices in Barnes et ai. (2015). 


Heller 


3. MODEL DESCRIPTION 

Our model evolves planet-star systems forward in time 
in order to determine whether HECs can form from mini- 
Neptunes that have migrated into the HZs of M dwarfs. 
We perform our calculations on a grid of varying plane¬ 
tary, orbital, and stellar properties in order to determine 
the types of systems that may harbor HECs. The com¬ 
plete list is provided in Table where we indicate the 
ranges of values we consider as well as the default values 
adopted in the plots in (unless otherwise indicated). 

Integrations are performed from t = tg (the time at 
which the planet is assumed to have migrated into the 
HZ) tot = tstop (the current age of the system) using the 
a daptive timestepping method described in Appendix E 
of [Barnes et al. ( [2013 ). 


3.1. Stellar Model 
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We use the evolutionary tracks of Baraffe et al. (1998 1 
for solar metallicity to calculate Lboi and 1 eff as a func- 
tion of time. We then use Q to calculate Exuv, given 
(l^xuv7^boi)sat = 10“^ and values of and j3 given in 
Table U 

Using Lboi and we calc ulate the location of th e 
HZ from the equations given in Kopparapu et al. (2013), 
adding the eccentricity correction H. Given the uncer¬ 
tainty in the actual HZ boundaries and their dependence 
on a host of properties of a planet’s climate, we choose 
our inner edge (IHZ) to be the average of the Recent 
Venus and the Runaway Greenhouse limits and our outer 
edge (OHZ) to be the average of the Maximum Green¬ 
house and the Early Mars limits. Throughout this paper 
we will also refer to the center of the HZ (GHZ), which 
we take to be the average of the IHZ and OHZ. Since 
we are concerned with the formation of ultimately hab¬ 
itable planets, we take the locations of the IHZ, GHZ, 
and OHZ to be their values at 1 Gyr, at which point the 
stellar luminosity becomes roughly constant. 

3.2. Planet Radius Model 

To determine the planetary radius Rp as a function 
of the core mass M^, the envelope mass fraction fn = 
Me/Mp, and the planet age, we use the p lanet struc¬ 
ture model descr ibed in Lopez et al. (2012) and |Lopez 
'fc Fortney (|2014), which is an extension of the model of 


h'ortney et ai. 12007) to low-mass low-density (LMLD) 
planets, 'i'hese models perform full thermal evolution 
calculations of the interior as a function of time. In our 
runs, the core is taken to be Earth-like, with a mixture 
of 2/3 silicate rock and 1/3 iron, and the envelope is 
modeled as a H/He adiabat. A grid of values of Rp 
is then computed in the range IM 0 < Me < 10 Mq, 
10"® < fn < 0.5, and lO^ears < t < 10^° years. For 
values between grid points, we perform a simple trilinear 
interpolation. For gas-rich planets, Rp is the 20 mbar 
radius; for gas-free planets, it corresponds to the surface 
radius. The evolution of Rp with age due solely to ther¬ 
mal contraction is plotted in Figure]^ for a few different 
core masses and values of fn- 
We note that the models of|Fortney et al.| (|2007|) and 
L opez et al. ( 2012 ) are in ge neral agreement wit h those 
otIMordasin i et ahl (|2012a|b|) and, by extension, |Rogers 
et al.|(|2011|). IMordakni et al.|(|201 2al) presented a valida - 
tion 01 their model against that of Fortney et al. (2007), 
showing that for planets spannihg 0.1 to 10 Jupiter 
masses, the two models predict the same radius to within 
a few percent. At the lower masses relevant to our study, 
the two models are also in agreement. To demonstrate 
this, in Figure we shade the regions corresponding to 
the spread in radii at a give n mass and age in Figure 
9 of Mordasini et al. (2012b). Since those authors used 
a coupled tormation/evolution code, at low planet mass 
the maximum envel ope mass fraction fn is s mall; for a 


total mass of 4 M 0 , Mordasini et al. (2012b) find that 


all planets have Jh < 0.2. At 2 M 0 , rnost planets have 
fn ^ 0.1. We can see from Figure that at these values 
of fn, the two models predict very sim ilar radius evolu¬ 
tion. Note that Mordasini et al. (2012b) did not consider 
planets less massive than 2 iVl 0 . 

The maximum envelope fraction merits further discus¬ 
sion. Since we do not model the formation of mini- 
Neptunes, we do not place a priori constraints on the 


.3 

cS 



<(yr) 


Fig. 3. — Evolution of the radius as a function of time due to 
thermal contraction of the envelope, in the absence of tidal effects 
and atmospheric mass loss. hVom top to bottom, the plots corre¬ 
spond to planets with initial total masses (core -h envelope) of 1, 
2, and 5 Mq. Line styles correspond to different initial hydrogen 
mass fractions: 1% (red, solid), 10% (green, dashed), 25% (blue, 
dot-dashed), and 50% (black, dotted). For comparison, the grey 
shaded reg ions in the bottom two p lots are the spread in radii cal¬ 
culated by|Mordasini et al.| (|2012b|| for fn ^ 0.20. See text for a 
discussion. 


value of Jh at a given mass; instead, we allow it to vary 
in the range 10"® < fn < 0.5 for all planet masses. At 
masses < 5 M 0 , planets accumulate gas slowly and are 
typically unable to accrete more than ~ 10 - 20 % of their 


mass in H/He 1 

Rogers et al. 

2011 

Bodenheimer & Lis-| 

sauer||2014 

); values of « 0 

.5 may thus be unphysical. 


around M dwarfs (Carpenter et al.||2006 Pascucci et al.j 


2009) allow more time tor gas accretion, potentially in- 
creasing the maximum value of fn- Nevertheless, and 
more importantly, if a planet with fjj = 0.5 loses its en¬ 
tire envelope via atmospheric escape, any planet with the 
same core mass and fn < 0.5 will also lose its envelope. 
Below, where we present integrations with fn = 0.5, 
our results are therefore conservative, as planets with 
Jh 0.5 will in general evaporate more quickly. 

While our treatment of the radius evolution is an im- 
provement upon pas t tidal-atmospheric coupling papers 
(Jackson et al.|2010 for instance, calculate Rp for super- 
Liarths by assuming a constant density as mass is lost), 
there are still issues with our approach: (1) We do not 
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account for inflation of the radius due to high insolation. 
Instead, we calculate our radii from grids corresponding 
to a planet receiving the same flux as Earth. While at 
late times this is justified, since planets in the HZ by 
definition receive fluxes similar to Earth, at early times 
this is probably a poor approximation; recall that plan¬ 
ets in the HZ around low mass M dwarfs are exposed 
to fluxes up to two orders of magnitude higher during 
the host star’s pre-main sequence phase. The primary 
effect of a higher insolation is to act as a blanket, de¬ 
laying the planet’s cooling and causing it to maintain an 
inflated radius for longer. This will result in mass loss 
rates higher than what we calculate here. (2) Since we 
are determining the radii from pre-computed grids, we 
also do not model the effect of tidal dissipation on the 
thermal evolution of the planet. Planets undergoing fast 
tidal evolution can dissipate large amounts of energy in 
their interiors, which should lead to significant heating 
and inflation of their radii. (3) The radius is also likely 
to depend on the mass loss rate. Setting Rp equal to the 
tabulated value for a given mass, age, and composition 
is valid only as long as the timescale on which the planet 
is able to cool is significantly shorter than the mass loss 
timescale. Otherwise, the radius will not have enough 
time to adjust to the rapid loss of mass and the planet 
will remain somew hat inflated, lead ing to a regime of 
runaway mass loss (Lopez et al.|2012). While the planets 


considered he r e are probably not in the runaway regime 


(Lopez et al. (2012) found that runaway mass loss oc- 
curred only tor H/He mass fractions > 90%), we might 


implications of this choice in §5.3[ 

Given the large planetary radii at early times, many of 
the planets we model here are not stable against Roche 
lobe overflow in the HZ. During RLO, the stellar gravity 
causes the upper layers of the atmosphere to suddenly 
become unbound from the planet; this occurs when > 
RRoche, where i?Roche is given by §. For a planet that 
forms and evolves in situ, RLO never occurs, since any 
gas that would be lifted from the planet in this fashion 
would have never accreted in the first place. However, 
an inflated gaseous planet that forms at a large distance 
from the star may initially be stable against overflow and 
enter RLO as it migrates inwards (since i?Roche cx a). 
This is particularly the case for planets in the HZs of M 
dwarfs, since a and consequently i?Roche are small. 

Ideally, the tidally-enhanced mass loss rate equation 
© should capture this process, but instead it predicts 
an infinite mass loss rate as Rxuv -RRoche (or as 
^ 1) and unphysically changes sign for ^ < 1. This 

is due to the fact that the energy-limited model implic¬ 
itly assumes that the bulk of the atmosphere is located 
at i?xuv (the single-layer assumption). Realistically, we 
would expect the planet to quickly lose any mass above 
the Roche lobe and then return to the stable hydrody¬ 
namic escape regime. However, upon loss of the material 
above i?Roche, the portion of the envelope below the new 
XUV absorption radius i?xuv ''^hl not be in hydrostatic 
equilibrium; instead, an outward flow will attempt to re¬ 
distribute mass to the evacuated region above, leading 
to further overflow. 

Several models exist that allow one to calculate the 


early active phase of the parent star. 

mass loss rate due to RLO (e.g., 

Ritter |1988 Trilling 

All points outlined above lead to an underestimate of 

et al.||1998 

Gu et al.||2003 

Sepinsky et al.| 

2U07 

). 'i'hese 


proportional to Rp (j^ or Rp'^^ (13), calculating the ra¬ 
dius in this fashion leads to a lower bound on the amount 
of mass lost and on the strength of the coupling to tidal 
effects. Because our present goal is to determine whether 
it is possible to form habitable evaporated cores via this 
mechanism, this conservative approach is sufficient. Fu¬ 
ture work will incorporate a self-consistent thermal struc¬ 
ture model to better address the radius evolution. 

3.3. Atmospheric Escape Model 

We assume that the escape of H/He from the planet at¬ 
mosphere is hydrodynamic (blow-off) at all time s, which 
is valid at t he XUV flu:^s we consider here (see Erkaev 
et al.||2013| and Figure [2|. We run two separate sets of 
integrations: one in which we assume the flow is energy- 
limited ([^ for all values of ExuV: one in which we 

switch fr om energy-limited to radiation/recombination- 
l imited ( |I^ above the critical value of the flux (see 
1 2.3.5|). fer planets whose orbits are eccentric enough 
that tney switch between the two regimes over the course 
of on e orbit, we make use of the expressions derived in 
the Appendix. These two sets of integrations 
oughly bracket the true mass loss rate. 

For eccentric orbits, we calculate the mass loss in the 
energ y-lim ited regime from (16), with ATecc evaluated 
from (1^). We vary exuv and Rxuv in the ranges given 
in Tal^lJ We choose exuv = 0.30 as our default case. 
While this is consistent with values cited in the literature 


change between the outflowing gas and the planet, which 
can lead to its outward migration, given by 


1 da 
a dt 


2 dMp 
Mr, dt 


( 22 ) 


for a planet on a circular orbit (Gu et al. 2003 Ghang 


for a 

planet 

et al. 

2010 ) 


(see 12.3.2), it could be an overestimate. We discuss the 


By differentiating the stability criterion Rx\jv{Mp) = 
i?Roche(A/p), one may then obtain an approximate ex¬ 
pression for dMp/dt in terms of the density profile dM{< 
R)/dR of the envelope. 

However, for mini-Neptunes that migrate into the HZ 
early on, RLO should occur during the initial migration 
process, which we do not model in this paper. Instead, 
we begin our calculations by assuming that our planets 
are stable to RLO in the HZ. If a planet’s radius initially 
exceeds the Roche lobe radius, we set its envelope mass 
equal to the maximum envelope mass for which it can be 
stable at its current orbit; the difference between the two 
envelope masses is the amount of H/He it must have lost 
prior to its arrival in the HZ. It is important to note that 
these planets will initially have Rxuv = RRoche, which as 
we mentioned above, leads to an infinite mass loss rate in 
(1^. An accurate determination of Mp in this case prob- 
Sly requires hydrodynamic simulations. However, the 
mass loss rate can be approximated by imposing a mini¬ 
mum value ^min in ([^ ■ For ^ , we set the mass loss 

rate equal to Mp{^ = ^min)- This is equivalent to impos- 
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ing a maximum mass loss enhancement factor l/il'tide, 
preventing the mass loss rate from reaching unphysical 
values as i?xuv —>■ -RRoche- 



10 

9 

8 

7 

6 

5 

4 

3 

2 

1 




Fig. 4. — Complete evaporation time tevap as a function of the 
cutoff value ^niin for a 2M0 mini-Neptune with fug = 0-5 on a- 
circular orbit around a 0.08 Mq star. The red, green, and blue lines 
correspond to planets in the IHZ, CHZ, and OHZ, respectively. 
Also plotted is the mass loss enhancement factor l/iftide (black 
dashed line), which approaches infinity as ^ —I- 1. Note that for 
?min iS 3, the evaporation time is relatively insensitive to the exact 
cutoff value, despite the fact that l/i^tide blows up. We therefore 
choose 5inin = 3 as the default cutoff, corresponding to a maximum 
enhancement factor l/lftide ~ 2. 


In Figure 1^ we show how the time tevap at which com¬ 
plete evaporation takes place scales with ^niin for a typ¬ 
ical mini-Neptune in the IHZ (red), CHZ (green), and 
OHZ (blue). Also plotted is the mass loss enhancement 
factor l/ATtide (Equation]^ black dashed line) as a func¬ 
tion of ^ = ^min- Interestingly, despite the fact that the 
instantaneous mass loss rate {Mp oc l/iCtide) approaches 
inhnity as ^ > 1, the evaporation time is relatively con¬ 
stant for ^niin < 3. This is because for very large Mp, the 
planet loses sufficient mass in an amount of time At to 
decrease Rp substantially and terminate the overflow. In 
other words, cases in which ^ « 1 are very unstable, and 
as mass is lost ^ will quickly increase beyond ~ 3. Both 
the net amount of mass lost and the evaporation time 
are therefore insensitive to the particular value of ^min) 
provided it is less than about 3. We therefore choose 
■^min = 3 as the default value for our runs, noting that 
this corresponds to a maximum mass loss enhancement 
of l/ATtide « 2. 


3.4. Tidal Model 

We calculate the evolution of the semi- majo r axis , the 
ecce n trici ty and the rotation rates from ( |B1[ )-(B3) and 
(Cl l-( C3) in the Appendix for the CPL and CTL models, 
respectively. For simplicity, we set the obliquities of all 
our planets to zero. Since the ti dal locking time scale 
is very short for close-in planet^ (|Gu et ah 2003), we 


^ Tidal locking refers to the state in which a body’s rotation rate 
is fixed by tidal forces at an equilibrium value. While for circular 
orbits, this implies uii = n, in the general case of an eccentric orbit 
in the CPL model, the plan et’s rotation rate as sumes a slightly 
super-synchronous value. See|Barnes et aL]||2013[l for a discussion. 


assume that the plan et’s rota tion rate is given by the 
equilibrium value ( |B6[ ) or ( |C6[ ). 

We calculate the stellar spin by assuming different ini¬ 
tial periods (see Table and evolving it according to the 
tidal equations, while enforcing conservation of angular 
momentum as the star contracts during the pre-main se- 
que nce phase. We neg lect the effects of rotational brak¬ 
ing (Skumanich 1972), whereby stars lose angular mo- 
mentum to winds and spin down over time. This leads 
to an overestimate of the spin rate at later times, but 
tidal effects should only be important early on, when 
the ra dii and the eccentricity are higher. [Bolmont et al.| 
(2012|) recently modeled the coupling between stellar spin 
and tides, following the evolution of a “slow rotator” star 
(Pq ~ 10 days) and a “fast rotator” star (Pg ~ 1 day). In 
both cases, the stars sped up during the first ^ 300 — 500 
Myr, after which time angular momentum loss became 
significant. However, tidal evolution is orders of mag¬ 
nitude weaker at such late times, so we would expect 
rotational braking to have a minimal effect on the tidal 
evolution. Moreover, as we show in the Appendix, in 
general it is the tide raised on the planet that dominates 
the evolution; as this depends on the planetary rotation 
rate, and not on the stellar rotation rate, our results are 
relatively insensitive to the details of the spin evolution 
of the star. 

In the CPL model, we adopt typical gas giant values 
10"^ < Qp < 10^ for gas-rich mini-Neptunes and typical 
terrestrial values 10 < Qp < 500 for planets that have 
lost their envelopes; we assume stellar values in the range 
10^ < Q* < 10®. In the CTL model, we consider time 
lags in the range 10“® s < Tp < 10^ s for gas-rich mini- 
Neptunes and 10~^ s < t„ < 10® s o nce t hey lose their 
envelopes. Following Leconte et al. (2010), we consider 
stellar time lags in the range 10“^ s < r* < 10“^ s. 

Once a mini-Neptune loses all of its atmosphere, we ar¬ 
tificially switch Qp or Tp to the terrestrial value adopted 
in that run. In reality, as the atmosphere is lost, the 
transition from high to low Qp (or low to high Tp) should 
be continuous. A detailed treatment of the dependence 
of Qp and Tp on the envelope mass fraction is deferred to 
future work. 

Finally, we note that the second-order CPL model 
described above is valid only at low eccentricity. For 
e > •\/l/19 « 0.23, the phase lag of the dominant tidal 
wave discontinuously changes from negative to positive, 
such that the model then predicts outward migration due 
to the planetary tide. This effect is unphysical, stemming 
from the fact that the CPL model considers only terms 
up to second order i n the eccentricity (fo r a detailed dis¬ 
cussion of this, see Leconte et al. 2010). We therefore 
restrict all our calculations m the CPL framework to val¬ 
ues of the eccentricity e < 0.2. For higher values of e, we 
use the higher-order CTL model. 


4. RESULTS 
4.1. A Typical Run 

In Figure]^ we show the time evolution of three mini- 
Neptunes as a guide to understanding the results pre¬ 
sented in the following sections. We plot planet mass and 
radius (top row), semi-major axis and eccentricity (cen¬ 
ter row), and stellar XUV luminosity and radius (bottom 
row), all as a function of time since the planet’s initial mi- 
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(a) Energy-Limited 


(b) Energy-Limited, to = 100 Myr (c) Radiation/Recombination Limited 
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Fig. 5.— Three sample integrations of our code. The first row plots show the envelope mass (left axis) and planet radius (right axis) 
versus time since formation; the second row plots show the semi-major axis (left) and eccentricity (right) versus time; and the third row 
shows the stellar XUV luminosity (left) and stellar radius (right) versus time. The planet is initially a 1 Mq core with a 1 M0 envelope 
orbiting around a 0.08 Mq M-dwarf with e = 0.5 at a semi-major axis of 0.04 AU, just outside the OHZ (light blue shading). Unless 
otherwise noted, all other parameters are set to their default values (Table [y. As it loses mass and tidally evolves, it migrates into the 
CHZ (light green shading). Note that the evolution of the HZ is not shown; the GHZ and OHZ are taken to be the long-term (> 1 Gyr) 
values, (a): In this run, we force the escape to be energy-limited, as in an X-ray dominated flow. The planet loses its entire envelope at 
t ^ 100 Myr. (b): Same as (a), except the calculation starts at to = 100 Myr, corresponding to a planet that undergoes late migration. 
While the envelope still completely evaporates, this occurs at a much later time, t ^ 2 Gyr. (c): Same as (a), except that the escape 
is radiation/recombination limited above the critical flux. Here, the envelope does not fully evaporate and tidal migration is noticeably 
weaker. See the text for a discussion of the labels A-F. 


gration, to, for 5 Gyr. In the center plots, the post-1 Gyr 
OHZ and GHZ are shaded blue and green, respectively. 
At t = to, the planet is a 1 M 0 core with a 1 Mq en¬ 
velope orbiting around a 0.08 Mq M-dwarf with e = 0.5 
at a semi-major axis of 0.04 AU. We set exuv = 0.30; 
other parameters are equal to the default values listed 
in Table [l] Because of the high eccentricity, the tidal 
evolution is calculated according to the GTL model. 

In column (a), we force the escape to be energy-limited 
([^, corresponding to an X-ray dominated flow. Prior 
to the first timestep, nearly 90 percent of the envelope 
mass is lost to RLO, indicated by the discontinuous drop 
marked A on the top plot. This is due primarily to the 
inflated radius shortly after formation, which reaches 30 
Rq for a planet of age t = tg = 10 Myr. Once this mass 
is removed, the planet enters the energy-limited escape 
regime, which operates on a timescale of ~ 10 Myr (B). 
After ~ 100 Myr (C), the planet loses its entire envelope 
and becomes a HEG. 

In the center plot, we see that the planet’s orbit 
steadily decays as it circularizes, with a sharp discontinu¬ 
ity in the slope at ^ 100 Myr (D), corresponding to when 
it transitions from a gaseous (low Tp) to a rocky (high 
Tp) body. The tidal evolution from that point forward is 
dramatically stronger, and e decreases to ^ 0.1 at 5 Gyr. 
The planet’s semi-major axis decays by 25%, moving it 
well into the GHZ. As we noted earlier, the transition 
from low to high tidal time lags (or, alternatively, from 
high to low tidal quality factors) is likely to be gradual 
as the bulk of the energy shifts from being dissipated in 
the envelope to being dissipated in the core. In this case, 


the faster inward migration as Tp increases is likely to ac¬ 
celerate the rate of mass loss, leading to slightly earlier 
evaporation times. However, given the large uncertainty 
in the values of Tp and its dependence on planetary and 
orbital parameters, our current approach should suffice. 

In the bottom plot, we see that the bulk of the mass 
loss occurs when the stellar XUV flux is high. After t « 
100 Myr, the XUV luminosity is low enough that a planet 
with significant hydrogen left {fn > 0.01) is unlikely to 
completely evaporate. Here, the XUV saturation time is 
set to 1 Gyr, visible in the kink marked by the label E; 
prior to that time, the decrease in the XUV luminosity 
is simply a function of the rate of contraction of the star. 
After t ^ 1 Gyr (F), the stellar radius asymptotes to 
its main sequence value and the XUV flux decays as a 
simple power law. 

In column (b) we repeat the integration but delay the 
start time, setting to = 100 Myr. This corresponds to a 
planet that undergoes a late scattering event, bringing it 
to a = 0.04 AU when both its radius and the XUV flux 
are significantly smaller. In this case, RLO is somewhat 
less effective, removing only 50% of the envelope initially 
(A). However, the planet still loses all of its hydrogen at 
t « 2 Gyr (C). Interestingly, because of its delayed evap¬ 
oration, the planet’s eccentricity at 5 Gyr is significantly 
higher than in the previous run. This occurs because the 
transition from low to high Tp (D) occurs much later. In 
this sense, a planet’s current orbital properties can yield 
useful information about its atmospheric history. How¬ 
ever, a more rigorous tidal model that accounts for the 
gradual change in Qp and Tp as fn decreases is probably 



























Habitable Evaporated Cores 


13 


necessary to accurately infer the atmospheric evolution 
based on a planet’s current eccentricity. 

Finally, in column (c) we repeat the integration 
performed in (a), but this time set the flow to be 
r adiatio n/recombination-limited above the critical flux 
(p.3.5|. Because of the significantly lower escape rate at 
early times, the envelope does not completely evaporate, 
and at 5 Gyr this is a super-Earth with slightly less than 
1% hydrogen by mass. In order for a planet to completely 
lose its envelope in the radiation/recombination-limited 
regime, it must either migrate into an orbit closer to the 
star, have a larger eccentricity, have a smaller core, or be 
stripped by another process. 



Initial (Mjg) Initial (Mjg) 


Fig. 6. — Initial versus final envelope mass (Mh) for planets 
that end up in the IHZ (red), CHZ (green), and OHZ (blue). 
Line styles correspond to different values of tg (solid, 10 Myr; 
dashed, 50 Myr; dash-dotted, 100 Myr). Columns correspond to 
runs in which the escape mechanism is energy-limited (left) and 
radiation/recombination-limited (right); rows vary certain param¬ 
eters as labeled, with all others set to their default values. In the 
default run, the planet has a 1 Mq core and orbits an M dwarf 
with M* = 0.08 Mq in a circular orbit. The dotted gray line corre¬ 
sponds to a planet that undergoes no evaporation. An X marks the 
critical initial envelope mass below which full evaporation occurs 
within 5 Gyr. In some plots, curves of a given color/line style are 
missing; for those runs, the entire envelope was lost for all starting 
values of Mjj. 


4.2. Dependence on Me, Mh, M*, and to 


In Figure we plot the initial versus final envelope 
mass [Mh) Iot planets that end up in the IHZ (red lines), 
CHZ (green lines), and OHZ (blue lines). Line styles cor¬ 
respond to different values of to, the time at which the 
planet migrated into its initial close-in orbit (solid, 10 
Myr; dashed, 50 Myr; dash-dotted, 100 Myr). The col¬ 
ored shadings are simply an aid to the eye, highlighting 
the spread due to different values of to- For reference, in 
dotted gray we plot the line corresponding to a planet 
that undergoes no evaporation. Note that in most plots, 
the curves approach this line as the initial Mh is in¬ 
creased: at constant core mass, it is in general more 
difficult to evaporate planets with more massive H/He 
envelopes. Finally, if a curve of a given color/line style 
is missing, the final hydrogen mass is zero for all values 
of the initial Mh- 

The two columns correspond to runs in which the 
escape was forced to be energy-limited (left) and 
radiation/recombination-limited at high XUV flux and 
energy-limited at low XUV flux (right). Rows correspond 
to different planet properties, as labeled. In the top (“De¬ 
fault”) row, the planet’s core mass is set to IM 0 . The 
planet orbits an M dwarf with M* = 0.08 Mq in a cir¬ 
cular orbit. All other parameters are set to their default 
values. 

In the second row, we double the core mass. The third 
row is the same as the first, but for a 0.16 Mq star; and 
the fourth row is the same as the first, but for an initial 
eccentricity e = 0.4. Since planets in this run undergo 
orbital evolution, the different color curves correspond 
to planets that end up in the IHZ, CHZ, and OHZ, re¬ 
spectively (their initial semi-major axes are somewhat 
larger). 

Let us first consider the general trends in the plots. 
For low enough Mh, all curves become nearly vertical. 
The critical envelope mass below which all mass is lost is 
marked with an X. Note that once planets lose sufficient 
mass such that Mh 10“"*, the envelope is very unsta¬ 
ble to complete erosion under the XUV fluxes considered 
here (corresponding to a near-vertical slope in this di¬ 
agram). Many curves also display a flattening towards 
high initial envelope masses; some have prominent kinks 
beyond which the final envelope mass is constant. This 
is due to Roche lobe overflow, which causes any planets 
with radii larger than the Roche radius to lose mass prior 
to entering the HZ. Since increasing the envelope mass 
increases the planet radius, this results in a maximum 
envelope mass for some planets. 

Planets that migrate in early (small to) lose signifi¬ 
cantly more mass than planets that migrate in late. This 
is a consequence of both the decay in the XUV luminosity 
of the star with time and the quick decrease of the planet 
radius as the planet cools. Another interesting trend is 
that the difference in the evaporated amount is much 
more pronounced in the energy-limited runs than in the 
radiation/recombination-limited runs. This is due to the 

Rp scaling of the mass loss rate in the latter regime 
(versus the steeper scaling in the former). The differ¬ 
ence in the initial radius across runs with different values 
of to is less significant in the radiation/recombination- 
limited regime, resulting in more comparable mass loss 

rates. We also note that the Fxjjv scaling of the mass 
loss rate in this regime results in a weaker dependence on 
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semi-major axis, as expected: the blue, green, and red 
curves (for a given to) are packed more closely together 
in the right column than in the left column. 

Now let us focus on individual rows. For the default 
run (a IM 0 core in a circular orbit around a 0.08 Mq 
star), all planets in the IHZ lose all of their hydrogen and 
form HECs, regardless of migration time, envelope mass, 
or escape mechanism. In the CHZ, only planets that 
migrate within < 50 Myr and undergo energy-limited 
escape form HECs for all initial values of Mh- However, 
HECs still form from planets with Mh ^ 0.5 — 0.9 Mq 
in the CHZ. In the OHZ, this is only possible for planets 
with less than about 1% H/He by mass. 

At twice the core mass (second row), all curves shift up 
and to the left, approaching the zero evaporation line for 
planets in the OHZ. In the IHZ, HECs still form from 
planets with any initial hydrogen amount for Iq = 10 
Myr and in the energy-limited regime. In all other cases, 
HECs only form from planets with Mh ^ O.I Mq. Of all 
the parameters we varied in our integrations, changing 
the core mass has the most dramatic effect on whether 
or not HECs can form. As we discuss below, habitable 
evaporated cores with masses greater than about 2 M 0 
are unlikely. 

At higher stellar mass (third row), HECs are 
again more difficult to form, particularly in the 
radiation/recombination-limited regime. Due to the 
more distant HZ, Roche lobe overflow is less effective 
in removing mass. The shorter super-luminous contrac¬ 
tion phase of earlier M dwarfs also results in less total 
XUV energy deposition in the envelope. However, in the 
energy-limited regime, HECs still form from planets with 
up to 50% H/He envelopes in the IHZ. 

Einally, the effect of a higher eccentricity (bottom row) 
is much more subtle. In general, these planets lose 
slightly less mass than in the default run, but the plots 
are qualitatively similar. At high eccentr icity, the orbit- 
integrated mass loss rate is higher (see (2.3.61, but be¬ 
cause of the orbital evolution, the planet must start out 
at larger a in order for it to end up in the HZ at 5 Gyr. 
These effects roughly cancel out: in general, habitable 
evaporated cores are just as likely on eccentric as on cir¬ 
cular orbits. Note, also, that the green curves in the bot¬ 
tom right plot are the only ones that are non-monotonic. 
The effect is very small, but hints at an interesting cou¬ 
pling between tides and mass loss. At high initial Mh, 
the radius is large enough to drive fast inward migra¬ 
tion (see below), exposing the planet to high XUV flux 
for slightly longer than a planet with smaller Mh (and 
therefore a smaller radius), resulting in a change in the 
slope of the curve at initial Mh ~ 0.15 M 0 . 


4.3. The Role of Tides 

Next, we consider in detail how tides affect the evolu¬ 
tion of HECs. We show in the Appendix that the net 
effect of tides is to induce inward migration and circular¬ 
ization of planet orbits in the HZ of M dwarfs, an effect 
that couples strongly to the atmospheric mass loss. Eor 
e < 0.7, the flux increases with time as planets tidally 
evolve, accelerating the rate of mass loss; at higher eccen¬ 
tricities, the flux actually decreases due to the circular¬ 
ization of the orbit (see the Appendix for a derivation). 
The changing mass and radius of the planet can then 
act back on the tidal evolution, either accelerating it (in 


cases where \dMp/dt\ ^ \dRp/dt\) or decelerating it in a 
negative feedback loop (otherwise). 

In Figurej^we show the results of an integration of our 
code on a grid of tidal time lag Tp versus the XUV ab¬ 
sorption efficiency exuv- Colors correspond to the final 
hydrogen mass fraction, f'pj] evaporated cores occur in 
the white regions {fn = 0). Dark gray indicates planets 
that either migrated beyond the HZ or remained exterior 
to it and are therefore not habitable. These plots pro¬ 
vide an intuitive sense of the relative importance of tidal 
evolution (j/-axis) and energy-limited escape (a;-axis) in 
determining whether or not a HEC is formed. We note 
that once a planet loses its gaseous envelope, we switch 
the time lag to r/ = max(rp, 64s), where the latter is 
a typical (ga s-free) tidal time la g of a rocky planet (see, 
for instance, Barnes et al.||20I^ . 

In the default run (aj, we consider a planet with a core 
mass of IM 0 and an envelope mass of IM 0 {fn = 0.5) 
orbiting at 0.07 AU in a highly eccentric (e = 0.8) or¬ 
bit around a 0.08 Mq star. The mass loss mechanism 
is radiation/recombination-limited escape at high XUV 
flux and energy-limited at low XUV flux. At a given 
value of log Tp, say 0 , the final hydrogen fraction is a 
strong decreasing function of exuVj as expected: the 
higher the evaporation efficiency, the smaller the final 
envelope. Interestingly, the dependence of fn on the 
tidal time lag can be nearly just as strong. At a given 
value of exuV) say 0.35, f'n depends strongly on r, vary¬ 
ing from 10“^ ® « 3% to 0%-that is, the tidal evolution 
controls whether or not a HEC forms. This is due to the 
fact that at large Tp, tidal migration is fast, bringing the 
planet into the HZ while its radius is still inflated and 
the stellar XUV luminosity is higher. Planets with lower 
values of Tp migrate in later and undergo slower mass 
loss. At very low Tp, tidal migration is too weak to bring 
the planets fully into the HZ. The effect is stronger for 
higher exuv^ these are planets whose radii decrease very 
quickly (due to the fast mass loss), slowing down the rate 
of tidal evolution and keeping them outside of the HZ for 
longer. 

In plot (b), we increase the core mass slightly to 
1.5 M 0 . The effect on the mass loss is significant, and 
HECs no longer form. At high exuv and high Tp, the 
lowest value of fn is about 1%. This reinforces what we 
argued above: HECs are most likely for the lowest mass 
cores. 

In plot (c), we instead decrease the eccentricity to 0.7. 
As in (b), HECs no longer form in these runs due to 
the decreased strength of the tidal evolution, and again 
the minimum fn is about 1%. Note that this does not 
mean that HECs are more likely at higher eccentricity in 
general—this is only the case here because we fix the ini¬ 
tial semi-major axis at 0.07 AU. At fixed final semi-major 
axis (i.e., what we can readily observe in actual systems), 
planets on initially circular orbits will still lose more mass 
than planets on initially eccentric orbits (which must 
have formed farther out). 



mass loss is significantly suppressed. 

Einally, in plot (e), we force the escape mechanism to 
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(b) Higher core mass 


(c) Lower eccentricity 


(a) Default run 
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Fig. 7.— Contours of the log of the hydrogen mass fraction f'^ at 5 Gyr as a function of the tidal time lag Tp and the XUV absorption 
efficiency exuv tor five integrations of our code. White regions correspond to planets that completely lost their envelopes; dark gray regions 
indicate planets that are not in the HZ at 5 Gyr. (a) The default run. The planet has a core mass of 1 Mq with /hq = 0.5, orbiting around 
a 0.08 Mq star in an initially highly eccentric orbit (e = 0.8) at a = 0.07 AU. All other parameters are set to their default values (TableQ, 
and the escape mechanism is radiation/recombination-limited at high XUV flux and energy-limited at low XUV flux. Note that the final 
envelope mass highly depends on both exuv ^rid Tp. (b) The same as (a), but for a higher core mass Me = 1.5Mq. No evaporated cores 
form in this scenario, (c) The same as (a), but for a lower eccentricity e = 0.7. Again, no habitable evaporated cores form, and the planet 
remains outside of the habitable zone for a larger range of Tp. (d) The same as (a), but for a shorter XUV saturation time of the parent 
star, fsat =0.1 Gyr. No evaporated cores form. Since the XUV flux drops off much more quickly, energy-limited escape is less effective in 
removing mass and thus is a weaker function of exuv- (®) The same as (a), but for energy-limited escape only, which would be the case 
if the flow is X-ray dominated. Note that in this case whether or not the planet becomes an evaporated core is a much stronger function 
of both Tp and exuv- 


be energy-limited only. Since the escape is now entirely 
controlled by exuv, the dependence on this parameter is 
naturally much stronger, and complete evaporation oc¬ 
curs in this case for exuv ^ 0.32 at any Tp. Because 
evaporation occurs more quickly in this case than in the 
other plots, at any given time the planet radii are smaller, 
resulting in less efficient migration at a given Tp. More 
planets therefore do not make it into the HZ. Interest¬ 
ingly, for very large exuv 0.4), all planets make it into 
the HZ. This is due to the fact that these planets tran¬ 
sition from gaseous (low Tp) to gas-free (high Tp, equal 
to 64s in these runs) very early on. Despite their lower 
radii, they benefit from the stronger tidal dissipation of 
fully rocky bodies and are able to make it into the HZ 
after 5 Gyr. 

In general, the coupling between tides and mass loss is 
quite complex. For planets with high initial eccentrici¬ 
ties, this coupling can ultimately determine whether or 
not a HEC will form. We discuss this in more detail in 



4.4. Evaporated Cores in the Habitable Zone 

Having shown that HECs are possible in certain cases, 
we now wish to explore where in the habitable zone we 
may expect to find them. For a “terrestrial” (we use 
this term rather looseljj^ planet detected in the HZ, it 
would be very informative to understand whether or not 
it may be the evaporated core of a gaseous planet, since 
its past atmospheric evolution may strongly affect its 
present ability to host life. To this end, we ran eight 
grids of 2.7 x 10® integrations each of our evolution code 
under different choices of parameters in Table For 
initial semi-major axes in the range 0.01 < oq < 0.5, ini¬ 
tial eccentricities in the range 0 < eo < 0.95, and stellar 
masses between 0.08 Mq < M* < 0.5 Mq, we calculate 

^ As we explain below, it is quite possible that HECs are all but 
terrestrial, since it is likely that they have large ice mass fractions 
and are compositionally distinct from Earth. By “terrestrial,” in 
this case, we mean planets that are similar in size and mass to Earth 
and are not gaseous; these may or may not have a rocky surface 
like Earth. Our definition of “terrestrial” therefore encompasses 
water worlds. 
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the final (i.e., at tstop = 5 Gyr in the default run) val¬ 
ues of a, e, and the envelope mass Mh- We then plot 
contours of Mh as a function of the stellar mass and 
the final semi-major axis and eccentricity (that is, the 
observable parameters). Since in principle, the mapping 
{aQ,eo, Mjjg) —>■ [a^e^Mn) is not necessarily a bijection 
(due to the nonlinear coupling between tides and mass 
loss and the finite resolution of our grid), the function 
Mij(a,e) may be multi-valued at some points. In these 
cases, we take Mh at coordinates (a, e) to be the mini¬ 
mum of the set of all final values of Mh that are possible 
at (o,e). Because of this choice, the Mh = 0 contour 
in a versus M* plots (Figures SpO) separates currently 
terrestrial planets that could be the evaporated cores of 
mini-Neptunes (left) from currently terrestrial planets 
that have always been terrestrial (right). In other words, 
we are showing where evaporated cores are ruled out. 

In Figure we perform the calculations described 
above for planets with Me = IM 0 and initial fn = 
0.5 (i.e., Mh = IM 0 ). The habitable zone is plot¬ 
ted in the background for reference, where again the 
IHZ, GHZ, and OHZ are indicated by the red, green, 
and blue shading, respectively. Line styles corre¬ 
spond to different final values of the eccentricity ( 0 , 
solid; 0.25, dashed; and 0.5, dash-dotted). Dark red 
lines correspond to the Mh = 0 contours in the 
radiation/recombination-limited escape model; dark blue 
lines correspond to the energy-limited model. Note 
that though we designate our two models as “energy- 
limited” and “radiation/recombination-limited”, at low 
XUV fluxes the escape is energy-limited in both models, 
since below Fcrit the energy-limited escape rate is always 
smaller than the radiation/recombination-limited escape 
rate. In this sense, the “radiation/recombination-limited 
model” is always conservative, as the mass loss rate is 
set to the minimum of Equations (§ and ( [T^ . 

As an example of how to interpret this figure, con¬ 
sider a IM 0 rocky planet discovered orbiting a 0.2 Mq 
star at O.I AU, that is, squarely within the GHZ. Since 
this planet lies to the left of all blue curves, under the as¬ 
sumption of energy-limited escape it could be a habitable 
evaporated core. On the other hand, if the atmospheric 
escape were radiation/recombination-limited, this planet 
must have always been terrestrial. 

Next, consider a putative rocky planet discovered 
around the same 0.2 Mq star skirting the inner edge 
of the HZ (i.e., at 0.07 AU). Under the energy- 
limited assumption, this planet could be a HEG. 
Eor radiation/recombination-limited escape, however, 
whether or not it could be a HEG depends on its present 
eccentricity. If the planet is currently on a circular orbit, 
we infer that it has always been terrestrial (since it lies 
to the right of the e = 0 contour). However, since the 
planet lies to the left of the higher eccentricity contours, 
if e > 0.25, it could be a HEG. 

We conclude from Figure that if energy-limited es¬ 
cape is the dominant mechanism around M dwarfs, plan¬ 
ets with Mp ~ IM 0 in the GHZ and IHZ of these 
stars could be habitable evaporated cores. For M* < 
0.15 Mq, HECs may exist throughout the entire HZ. If, 
on the other hand, these planets are shaped mostly by 
radiation/recombination-limited escape, HEGs are only 
possible in the IHZ for M* < 0-2 and in the GHZ of 
M dwarfs near the hydrogen-burning limit. 


The effect of the eccentricity is significant primarily 
in the radiation/recombination-limited regime, where a 
lower mass loss rate keeps the planet’s radius large for 
longer than in the energy-limited regime, allowing it to 
migrate into the HZ faster, thereby enhancing the mass 
loss rate. In some cases, particularly near the inner edge 
of the HZ, the present-day eccentricity yields important 
information about a planet’s evolutionary history: de¬ 
pending on the precise value of e, a given terrestrial 
planet may or may not be a HEG. However, we urge 
caution in interpreting the results in Figure!^ at nonzero 
eccentricity, given the large uncertainty in the tidal pa¬ 
rameters of exoplanets. 

On this note, it is important to bear in mind that the 
curves in Figure!^ are a function of our choice of parame¬ 
ters in Table[l] order to assess the impact of our choice 
of “default” parameters on these contours, in Figure 
we repeat the calculation for more conservative values 
of the two parameters that govern the mass loss mecha¬ 
nism: the XUV saturation time fgat and the absorption 
efficiency exuv- In this figure, we choose Gat = 0-1 Gyr, 
the nominal value for earlier K/G dwarf stars ((2.1), and 
exuv = 0.15. 

In this grid, all curves shift quite dramatically to 
lower a, and HEGs are no longer possible under 
radiation/recombination-limited escape. Eor energy- 
limited escape, HEGs are confined to the IHZ for M* 
0.15 Mq and the GHZ for AG G O.I Mq at all eccentrici¬ 
ties considered here. 

Given the large difference between the results of Fig¬ 
ures 1^ and [9 care must be taken in assessing whether a 
planet ma 5 rbe a HEG. Since it is likely that the XUV 
saturation time is much longer for M dwarfs than for K 
and G dwarfs, and since our goal at this point is to sep¬ 
arate regions of parameter space where HECs can and 
cannot exist. Figure!^ is probably the more relevant of 
the two. We discuss this point further in ^ 

In Figure [T0| we raise the core mass toMc = 2 M 0 . 
HECs are now possible only in the IHZ and only in the 
energy-limited regime. At the lowest M*, HECs may 
be possible in the GHZ (energy-limited) and at the very 
inner edge of the HZ (radiation/recombination-limited). 
At higher eccentricity, the parameter space accessible to 
HECs is slightly higher in the energy-limited regime, but 
it is still a small effect overall. Eor a conservative choice 
of escape parameters (Gat = O.I Myr and exuv = 0.15) 
with Me = 2 M 0 , no HECs form anywhere in the HZ. 

5. DISCUSSION 
5.1. Initial Conditions 

We have shown that HECs can form from small mini- 
Neptunes {Mp < 2 M 0 ) with H/He mass fractions as high 
as fn = 0-5 around mid- to late M dwarfs. However, 
instrumental to our conclusions is the assumption that 
these planets formed beyond the snow line and migrated 
quickly into the HZ at t = G- Planets that form in situ 
in the HZ are unlikely to accumulate substantial H/He 
envelopes due to large disk temperatures and relatively 
long formation timescales. Significant gas accretion can 
only occur prior to the dissipation of the gas disk, which 
occurs on a timescale of a few to ~ 10 Myr; in particular. 


hammer et al. (20141 showed that a 1 M 0 terrestrial 
planet generally accretes less than ~ 3% of its mass in 
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Default Run, 1 M 0 



Fig. 8.— Regions of parameter space that may be populated by HECs, for Me = IM 0 , Mh < IM 0 , and default values for all other 
parameters. Terrestrial planets detected today occupying the space to the left of each contour line could be the evaporated cores of gaseous 
planets with fu < 0.5. Planets detected to the right of the contour lines have always been terrestrial/gaseous. Dark red lines correspond to 
the conservative mass loss scenario, in which mass loss is radiation/recombination-limited at high XUV flux and energy-limited at low XUV 
flux. Dark blue lines correspond to mass loss via the energy-limited mechanism only. Planets around stars with significant X-ray emission 
early on are likely to be in the latter regime. Different line styles correspond to different eccentricities today. Terrestrial planets detected 
at higher eccentricity (dashed and dash-dotted lines) could be evaporated cores at slightly larger orbital separations than planets detected 
on circular orbits (solid lines). Note that in the energy-limited regime, all 1 Mq terrestrial planets in the HZ of low-mass M dwarfs could 
be habitable evaporated cores. At higher stellar mass, HECs are restricted to planets in the CHZ and IHZ. In the radiation/recombination- 
limited regime, the accessible region of parameter space is smaller, but around the lowest mass M dwarfs HECs are still possible in the 
CHZ. 


H/He in the HZ of a solar-type star. Moreover, planets 
that form in situ do not undergo Roche lobe overflow, 
since accretion of any gas that would lead to Rp > i?Roche 
simply will not take place. 

But can mini-Neptunes easily migrate into the HZ? 
The large number of recen tly discovered hot Neptunes 
and hot Super-Earths (e.g., Howard et al. ]M 2 | ) suggests 
that inward disk-driven migration is an ubiquitous pro¬ 
cess in planetary systems, a s it is highly unlikely tha t 
these systems formed in situ (|Raymond & Cossouj2014|). 
Moreover, systems such as GJ 18U, GJ 422, and GJ 667C 
each have at least one super-Earth/mini-Neptune in the 
HZ (Tuomi et al. 2014 Anglada-Escude et al. 12013 ), some 
or ail ot which may haive migrated to their current orbits. 
On the theor etical front, N-body simulations by |Ogihara| 
& Ida (2009|) show that migration of protoplanets into 
the HZ ot M dwarfs is efficient due to the proximity of 
the ice line and the fact that the inner edge of the disk 
lies close to the HZ. Mini-Neptunes that assemble early 
beyond the snow line could in principle also migrate into 
the HZ, but the migration probability and the distri¬ 
bution of these planets throughout the HZ needs to be 
investigated further in order to constrain the likelihood 
of formation of HEGs. 

A second issue is whether or not Earth-mass cores can 
accrete large H/He envelopes in the first place. One of 
our key results is that Earth-mass HEGs can form from 
mini-Neptunes with initial envelope mass fractions of up 
to 50%; this would require a 1 M® planetary embryo to 
accrete an equivalent amount of gas from the disk. While 


such a planet would be stable against RLO beyond the 
snow line, whether or not it could have formed is un¬ 
clear. Gas accretion takes place in two different regimes: 
(i) a slow accretion regime, in which the envelope re¬ 
mains in hydrostatic equilibrium and gains mass only as 
it cools and contracts, evacuating a region that is then 
filled by nebular gas; and (ii) a runaway accretion regime, 
in which the rapidly increasing mass of the envelope leads 
to an increase in the size of t he Hill sphere and i ncreas- 
ingly faster gas accretion (e.g., [Pollack et al.|1996 ). Since 
the critical core mass for runaway accretion is thought to 
be somewhere in the vicinity of 10-20 M 0 (Pollack et al. 


1996 Rafikov 2006), the progenitors of HEGs must ac¬ 
crete 0 eir g as slowly. As we mentioned above, ILammer] 


et al. (2014) find that Earth-mass planets typically form 
with jH ^ 0.03 in the HZ of solar-type stars. This is 
consistent wit h rec ent analyses of Kepl e r plan et data by 


Rogers (2014) and Wolfgang & Lopez (2014), who find 
that most planets with radii less than about 1.5 R 0 (cor¬ 
responding to masses less than ~ 5 M 0 ) are rocky, with 
typical H/He mass fractions of about 1%. 

Eormation beyon d the snow line can increase fn ■ 


Rogers et al. (2011) show that a planet with core mass 
2.65 MfB accret es only 0.54 Mm of gas ( corres ponding 
to fn ~ 17%). |Bodenheimer & Lissauer| (|2014|), on the 
other hand, show that planets with cores m the range 
2 . 2 - 2 .5 M 0 generally accrete less than 10 % of their mass 
in H/He. However, these authors terminated gas a ccre¬ 
tion at 2 Myr. Eor a l o nger d isk lifetime of 4 Myr, Bo- 
denheimer & Lissauer (2014) show that a planet with 
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Conservative, 1 M0 



Fig. 9. — The same plot as Figurel^ but for a conservative choice of the parameters governing atmospheric escape: a short XUV saturation 
time tsat = 0-1 Gyr and a low XUV absorption efficiency exuv = 0.15. In this case, HECs are no longer possible for radiation/recombination- 
limited escape. For energy-limited escape, HECs are only possible in the IHZ of low-mass M dwarfs and in the CHZ of M dwarfs at the 
hydrogen-burning limit. At high eccentricity (e = 0.5), HECs are only marginally more likely. 


core and envelope masses of 2.8 M® and ~ 1.2 M^, re¬ 
spectively, can form {fn ~ 30%). This is particularly 
relevant to planet formation around M dwarfs, as these 


ter et al. 

2006 

tor larger 

initial 


Pascucci et al. 2009), which could allow 
H/He envelope tract ions. 

Nevertheless, while a 1 M 0 core with a 1 M 0 envelope 
is probably an unlikely initial condition, it is important 
to keep in mind that our results apply to planets that 
form with smaller H/He fractions as well. Figures 
10 show where planets with up to 50% H/He can form 
HECs; planets with lower H/He mass fractions also form 
HECs at the same values of a and M*. See §5.6| for a 
more detailed discussion. 

5.2. Are HECs Habitable? 

Under the core accretion mechanism, mini-Neptunes 
are likely to form close to or beyond the snow line, where 
disk densities are higher due to the condensation of var¬ 
ious types of ices; these planets should therefore have 
large quantities of volatiles, including water, ammonia, 
and CO 2 ices. If we assume a disk composition similar 
to that around the young Sun, the bulk composition of 
their cores would probably be similar to that of comets: 
roughly equal parts ice and silicate rock, similar t o what 
studies predict for the compos ition of water worlds (Leger 
et al. 2004 Selsis et al. 2007). Once these planets have 
migrated into the HZ and lost their envelopes, it is quite 
likely that they would be water worlds and therefore not 
“terrestrial” in the strictest sense of the word. Whether 
or not such planets are habitable may depend on their 
ability to sustain active geochemical cycles, which are 
crucial for life on Earth today. 

One concern is a possible interruption of the carbon 
cycle by a high pressure ice layer at the bottom of the 


ocean (Leger et al. 2004). Recently, Alibert (2014) cal¬ 
culated tEeEHtlcar ocean mass for high pressure ice for¬ 
mation, finding that it lies between 0.02 and 0.03 M 0 
for an Earth-mass planet. If HECs are in fact comet¬ 
like in composition, a deep ice layer would separate their 
oceans from their mantles, which could inhibit the recy¬ 
cling of carbon and other bioessential elements between 
the two reservoirs, a process that is critical to life on 
Earth. However, whether this is the case is far from 
settled, as processes involving solid state ice convection 
could mediate th e cycli ng of these elements. In partic- 


ular, Levi et al. (2013) and Kaltenegger et al. (12013) 


presented a mechanism that could recycle CH 4 and CU 2 
between the interior and the atmosphere of water worlds, 
invoking the ability of these molecules to form clathrates 
that could be convectively transported through the ice 
to the surface. Other mechanisms could also exist, and 
without further modeling, it is unclear whether a high 
pressure ice layer poses a threat to habitability. 

The probable difference in composition between HECs 
and Earth is likely to have other geophysical implica¬ 
tions. Hydrogen-rich compounds such as methane and 
ammonia could make up a non-negligible fraction of a 
HEC’s mass, leading to extremely reducing conditions at 
the surface. These could also end up in a secondary at¬ 
mosphere along with large quantities of CO 2 , which could 
lead to strong greenhouse heating, although it is likely 
that most of th e CO 2 and NHcj wo uld be sequestered in 
the ice mantle (Leger et al. 2004). The compositional 
difference of HECs would also likely lead to mantle con¬ 
vection and tectonic activity different from Earth’s, as 
well as differences in magnetic field generation by a pos¬ 
sible core dynamo. Since both an active tectonic cycle 
and a magnetic field may be necessary for habitability, 
these issues need to be investigated further. 
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Default Run, 2 M 0 



Fig. 10.— The same as Figure]^ but for a higher core mass Me. = 2 Mq. HECs are now confined mostly to the IHZ of low-mass M dwarfs 
in the energy-limited regime. For planets undergoing radiation/recombination-limited escape, super-Earth HECs may not be possible. 


Also critical to the habitability of a HEC is its abil¬ 
ity to outgas a secondary atmosphere once its primordial 
envelope is lost. In particular, the secondary atmosphere 
must be stable against erosion. While the power-law de¬ 
cay in XUV emission after fgat could allow for such an at¬ 
mosphere to form, low mass M dwarfs may remain active 
for t 1 Gyr. Strong planetary m agnetic fields could 
potentially shield these atmospheres; Segura et al. (20101 
showed that flares on the extremely active M dwarf AD 
Leo would have a small effect on the atmosphere of an 
Earth-like planet in the HZ provided it has an Earth¬ 
like magnetic field. However, interactions with the stel¬ 
lar wind could still pose serious problems to these and 
other planets in the HZs of M dwarfs. In particular, 
the high XUV/EUV fluxes of active M dwarfs can lead 
to significant expansion of the upper atmosphere, poten¬ 
tially causing the radius of the exo base to exceed the 
distance to the stellar mag netopause (Lichtenegger et al. 
20101 Lammer et a l. 2011 ). For a pure IN 2 atmosphere, 
Lichtenegger et al.| (|201Up showed that nitrogen ionized 
by EUV radiation above the exobase is subject to ion 
pickup by the solar wind, leading to the complete erosion 
of a 1 bar atmosphere in as short a time as 10 Myr. A 
stronger planetary magnetic field or large quantities of a 
heavier background gas such as CO 2 may be necessary to 
suppress ion pickup and preserve secondary atmospheres 
on HECs. 

Given the complex processes governing the habitability 
of HECs, a detailed investigation of these issues is left to 
future work. As we discuss below, stronger constraints 
on the details of the X-ray/EUV evolution of M dwarfs 
of all masses are essential to understanding the fate and 
ultimately the habitability of planets around these stars. 

5.3. The Need For Better Constraints 

Figure shows that the formation of HECs de¬ 
pends strongly on the atmospheric escape mechanism. 


As we mentioned earlier, whether a flow is closer to 
radiation/recombination-limited or energy-limited will 
depend on the ratio of the X-ray to EUV luminosity of 
the parent star. Low-mass M dwarfs may have XUV lu¬ 
minosities as high as ~ 4 X 10^® erg/s early on (^2.1). If 
X-rays contribute more than a few percent of this fumi- 
nosity, low-mass low-density planets in the HZs of these 
stars may un dergo energy-limited X- ray-driven escape 
(Figure 11 in Owen & Jackson 2012). Unfortunately, 
knowledge of the exact age-fumihosity relation in the X- 
ray and EUV bands for M dwarfs is still very poor, in 
great part because of the large uncertainties on these 
stars’ ages. However, recent studies suggest that X-rays 
contribute a significant fraction of this luminos ity (for 


Scalp et al. 2007). In particular, Stelzer] 
et al. (2013) report high {Lx > 10^® erg/s) X-ray lu- 


a review, see 


mmosities for a sample of early active M dwarfs in the 
solar neighborhood, with a steeper age dependence than 
in EUV and NUV bands, which dominate the emission 


for t > I Gyr. This is consistent with Owen & Jack- 
(2012), who argue that close-in mini-JNeptunes may 


son 


undergo a transition from X-ray-driven escape at early 
times to EUV-driven escape at later times. 

Moreover, atmospheric X-r ay heating and cooling is 
primarily done by metals. As Owen & Jackson (2012) 
point out, atmospheric composition is also likely to play 
a role in determining whether hydrodynamic flows are 
EUV- or X-ray-driven. Additionally, the presence of dust 
in the envelopes of these planets could greatly affect their 
ability to cool via Lyman a radiation. Absorption of re¬ 
combination radiation by dust particles lifted high into 
the envelope by vigorous convection could convert a sig¬ 
nificant fraction of this energy into heating, which could 
effectively increase the absorption efficiency exuv and 
bring the flow closer to the energy-limited regime. Un¬ 
fortunately, our present parametric escape model is un¬ 
able to address the effect of dust and metal abundances 
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on the escape rate—this issue needs to be revisited in the 
future with more sophisticated hydrodynamic models. 

The difference between Figurej^and Figurej^is just as 
signihcant. At lower tgat and exuV; the H EC p arameter 
space is greatly reduced. As we discuss in §2.11 the lower 
value tgat =0.1 Gyr is more representative of K and G 
dwarfs, while M dwarfs may remain saturated for tgat > 1 
Gyr. In fact, the energy-limited contours in Figure 
closely trace the CHZ/OHZ boundary as M* increases, 
predicting that HECs may even be possible around solar 
type stars (not shown). This is of course not the case, 
since solar-type stars are kn own to leave the saturation 
phase around Gat = 0.1 Gyr (Ribas et al. 2005). For some 
stellar mass near the M/K Dwarf transition (0.5 Mq < 
^ 0.7 Mq), Figure [^becomes a better description 
for HECs, but for low-massTvI dwarfs Figurej^is probably 
more appropriate. 

Naturally, the exact value of exuv '''^ill also affect the 
possible distribution of H ECs within the HZ. Recently, 


Shematovich et al. (2014) performed numerical simula- 
tions to solve the kinetic Boltzmann equation for XUV- 
irradiated hydrogen atmospheres, finding an upper limit 
to the heating efficiency of exuv ~ 0.20. Our nominal 
value exuv = 0.3 may therefore be an overestimate for 
many mini-Neptunes. However, given that our goal in 
this paper is to explore where in the HZ HECs are possi¬ 
ble and to map regions where the transition from gaseous 
to terrestrial is not possible, our present approach should 
suffice. Nevertheless, further studies constraining exuv 
are crucial to improving our understanding of this pa¬ 
rameter space. 


5.4. Eccentricity Effects 

Erom Figures |8 |[To 1 we see that as we consider planets 
with higher current eccentricity, the region where evapo¬ 
rated cores are possible overlaps increasingly more of the 
HZ, particularly in the radiation/recombination-limited 
regime. This does not necessarily mean that HECs are 
more likely at high eccentricity—but rather that plan¬ 
ets in circular orbits at certain a cannot be HECs, while 
planets at higher eccentricity can. However, that being 
said, there is an interesting negative feedback between 
tides and mass loss that may enhance the probability of 
HECs on eccentric orbits. A planet with initially high 
eccentricity will undergo fast tida l ev olution (due to the 
strong p dependence in Equation |C2[ ), particularly if its 
radius is large (since de/dt oc Rp). If a planet loses mass 
quickly (as is usually the case for gaseous planets with 
large Rp and at large eccentricity), its radius will shrink 
and de/dt will decrease. The earlier a planet sheds its 
envelope, the more likely it is to maintain a nonzero ec¬ 
centricity in the long run (i.e., after ~ 5 Gyr). Since 
HECs typically form from such quickly evaporating plan¬ 
ets, they are more likely to end up in eccentric orbits 
than their gaseous counterparts that did not lose their 
envelopes. There is, of course, a trade-off here in the 
sense that gas-free planets should have higher Tp (and 
therefore faster de/dt) than gaseous planets, but the de¬ 
pendence on Tp is linear and therefore much weaker than 
the dependence on the radius. What this means is that 
there is an interesting link between present eccentricity 
and mass loss history. Translating a planet’s present ec¬ 
centricity into a probability that it is an evaporated core 


is no easy task, however, and likely requires a detailed un¬ 
derstanding of its initial orbital state and the migration 
mechanism(s) it underwent. On the other hand, statis¬ 
tical surveys of planets found to the left of the contours 
in Eigure|^ could uncover interesting trends. 

We note that Eigures [8p0| show eccentricity contours 
only up to e = 0.5. At fmal eccentricities higher than 
0.5 today, the contour lines should move farther to 
the right, increasing the parameter space populated by 
HECs. However, since high eccentricities today in gen¬ 
eral require extremely high eccentricities in the past, we 
do not calculate contours above this level. We note, fi¬ 
nally, that even thoug h w e run simulations with eg as 
high as 0.95, Eigures SpO remain unchanged for a lower 
cutoff Co < 0.7. 


5.5. Orbital Effects Due To Roche Lobe Overflow 

In the present work we neglect any orbital effects due 
to Roche lobe overflow, which could in some cases lead 
to significant outward migr ation of the planet; see (p 2 ). 
We argued in section 


_I that since RLO should oc¬ 
cur during the initial stage of migration (i.e., from be¬ 
yond the snow line into the HZ), modeling it was outside 
the scope of the paper. However, a careful investigation 
of this initial migration process could uncover interest¬ 
ing couplings. Eor instance, mini-Neptunes that initially 
overshoot the HZ would experience strong RLO, which 
could in principle cause them to migrate back into the 
HZ. These planets could have lost significantly more mass 
than the ones we considered here, since both RLO and 
XUV-driven escape would be stronger in their closer-in 
orbits. 

A second important point concerning RLO -induced mi¬ 
gration is that for nonzero eccentricity, ( 22 ) does not ap¬ 
ply. In this case, it can be shown that angular momentum 
exchange between the gas and the planet at pericenter 
(where overflow is strongest) leads to a net increase in 
the eccentricity. In the limiting case that both bodies 
may be treated as point masses and the mass transfer 
rate may be approximated as a delta f unction at peri¬ 
center (which is appropriate for high e) 

(2009| show that 


Sepinsky et al. 
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(23) 

(24) 


predicting that both a and e will tend to increase with 
time (recall that Mp is negative). Their relative rates of 
change are 


de ,1 da 

“u = (1 ~ e)——. 

dt a dt 


(25) 


In other words, at intermediate values of e, the fractional 
rates of change in the semi-major axis and the eccen¬ 
tricity are comparable, and the eccentricity will increase 
proportionally to the semi-major axis. 

As e increases, so does the atmospheric mass loss rate 
(via ATecc), the Roche lobe overflow rate (modulo the in¬ 
crease in a), and the rate of tidal evolution, leading to 
potentially faster mass loss and rich couplings between 
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these different processes. We will consider all these ef¬ 
fects in a future paper. 

5.6. Different fn 

Figures [8p0| correspond to planets with initial hydro¬ 
gen fractioii fn = 0.5. Planets with lower initial fn 
are naturally more unstable to complete loss of their 
envelopes—this would move all contours to the right, 
expanding the region where HECs are possible. How¬ 
ever, for a different choice of initial fn in the range 
0-1 ^ Jh ^ 0.5, the figures change very little: at constant 
core mass, it is only marginally harder to fully evaporate 
du fn = 0.5 envelope than it is to evaporate & Jh = 0.1 en¬ 
velope (see discussion below). This is due to both Roche 
lobe overflow, which dramatically reduces the envelope 
mass early on, and the fast atmo^heric escape for ex¬ 
tremely inflated planets. In Figurel^a), for instance, the 
envelope fraction decreases by a factor of 10 within the 
first ^ 10 Myr due to energy-limited escape. The escape 
process is in general very fast for large fm and bottle¬ 
necks for fn fS 0.1; thus the HEC boundary is relatively 
insensitive to the exact choice of the initial fn- 



Fig. 11.— Similar to Figure]^ but here contours correspond to 
different choices of final fu for a 1 M 0 core on a circular orbit. The 
solid lines correspond to f'„ = 0 and are the same as the e = 0 
contours in Figure Dashed lines correspond to the transition 
between planets that have less than (left) and more than (right) 
0.1% H/He = 0.001) at 5 Gyr. Dotted lines correspond to the 
1% H/He (f'jj = 0.01) transition. While the f'jj = 0.001 contours 
are barely distinguishable from the = 0 contours, at final = 
0.01 the HEC parameter space is significantly larger. However, it 
is unclear whether planets with f'^ = 0.01 could be habitable. 


Along the same lines, we can ask how the choice of 
final Jh affects our results. In this study, we define 
an evaporated core as a planet with = 0 and no 
atmosphere. However, mini-Neptunes with substantial 
gaseous envelopes that evaporate down to f'^j ~ 0.01 be¬ 
come fundamentally solid planets. At such low /^, addi¬ 
tional non-thermal mass loss processes, including flaring 
events and interactions with the stellar wind, may signif¬ 
icantly erode the remaining envelope. Investigating the 
habitability of planets with nonzero is beyond the 
scope of this paper, but it is worthwhile to consider how 


an alternative definition of a “habitable evaporated core” 
affects our results. To this end, in Figure [^we plot con¬ 
tour lines corresponding to three different choices of the 
critical below which we consider a planet to be a HEC. 
Solid lines correspond to /^ = 0 (the default choice); 
dashed lines correspond to = 0.001; dotted lines, to 
f'jj = 0.01. Note that in the radiation/recombination- 
limited regime, this choice does not significantly affect 
evaporated cores in the HZ. However, for = 0.01, 
all planets undergoing energy-limited escape in the HZ 
could be HECs. 

Earths and super-Earths with up to a few percent of 
their mass in H/He may still harbor liquid water oceans, 
and so at this point their habitability can not be ruled 
out. At higher //f, how ever, the surface pressures in 
excess of ~ 1 GPa (e.g., Choukroun & Grasset 20071 
result in the formation of high-pressure ices, at which 
point the planet will likely no longer be habitable. 

5.7. Other Atmospheric Escape Mechanisms 

We modeled the XUV emission of M dwarfs as a 
smooth power law, but active M dwarfs are seldom so 
well behaved. Frequent flares and coronal mass ejections 
punctuate the background XUV emission and will lead 
to erosion of planet atmospheres beyond what we model 
here. During flares, the ratio of the X-ray to bolometric 
luminosity (L x /Trot,) can incre ase by up to two orders 
of magnitude ( Scalo et al.||2007 ). Moreover, interactions 
with the stellar wind and nonthermal escape mechanisms 
such as ion pick-up and charge-exchange can lead to a few 
Earth-ocean equivalents of hy drogen from t hese planets 
dKislyakova et ah 2013, 20141. However, as Kislyakova 
et al. 12013) demonstrate, the escape rate due to these 


processes generally amounts to only a few percent of the 
hydrodynamic escape rate. Including these non-thermal 
mechanisms would thus have a minor effect on our re¬ 
sults. 

In this study we also neglect the effects of magnetic 
fields. The presence of a strong planetary magnetic field 
may inhibit escape during flares and possibly decelerate 
the mass loss rate if the planetary wind is ionized. Since 
a strong magnetic field is likely a requirement of su rface 
habitability around an M dwarf ( Scalo et al.||2007 ), this 
point must be addressed in future work. 

Finally, stars less massive than about O.IMq may be in 
a “supersaturation” regime early on, saturating at XUV 
fluxes one or eve n two orders of ma gnitude below higher 
mass M dwarfs (Cook et al. 2014). This could reduce 
escape rates froni planets around M dwarfs close to the 
hydrogen burning limit. 

5.8. Multi-Planet Systems 

In this paper we restricted our calculations to single¬ 
planet systems. However, many of the super-Earths 
and mini-Neptunes recently detected by Kepler reside in 
multi -planet syste ms, such as Kepler-3 2 (Fabrycky et al. 
2 OI 2 I), Kepler-62 (iBoru cki et al. |20 13|), arid Kep l er-186 
iGuintana et al. 2014 Moreover, Swift et al. (2013) 
demonstrated how the live-planet system Kepler-32 could 
be representative of the entire ensemble of Kepler M 
dwarf systems, arguing that multiplicity could be the 
rule rather than the exception for these stars. 

While tidal processes still operate in systems with mul¬ 
tiple planets, the orbital evolution of individual planets 
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will be much more complex. In general, though tidal 
dissipation may still lead to a decrease in the semi¬ 
major axis, the eccentricity evolution will be governed 
by secular interactions between the planets. The cou¬ 
pling betw een mass loss and tidal evolution we inves¬ 
tigated in §4.3| would likely be weaker, particularly for 
closely-packed coplanar systems where the eccentricities 
are necessarily low. However, planet-planet interactions 
could allow for even richer couplings to the atmospheric 
evolution of planets in the HZ. 

Consider, for instance, the case of a terrestrial planet 
in the HZ and a mini-Neptune interior to the HZ. Both 
mass loss and tidal dissipation scale inversely with the 
semi-major axis, so both effects could strongly shape the 
fate of the inner planet, which could in turn affect both 
the orbital and atmospheric evolution of the HZ planet 
in complex ways. The same could happen for a planet 
overflowing it s Ro che lobe interior to the HZ, which as 
we argued in §5.5| could experience large changes in a and 
e. 


Similarly, consider a scenario in which a planet in the 
HZ loses mass at the same time as it is perturbed by 
a more massive companion exterior to the HZ. As the 
planet’s mass decreases, it will undergo larger swings in 
eccentricity, which in turn lead to faster mass loss and 
stronger orbital evolution. 

Another interesting scenario involves planets close to 
mean moti on resonances. F o r a sy stem of two planets in 


resonance. 


Lithwick & Wu| (f20l^ and |Batygin & Mor- 
:ne 


bidelli (20l3p showed that the presence of a dissipative 


mechanism such as tidal evolution naturally leads to the 
repulsion of the two planets’ orbits: the inner planet’s or¬ 
bit decays, while the outer planet gets pushed outward. 
In the hypothetical case of two mini-Neptunes in res¬ 
onant orbits interior to the HZ, this mechanism could 
result in the migration of the outer planet into the HZ, 
having lost significantly more mass (due to its initially 
smaller orbit) than if it had originated exterior to the 
HZ and migrated inwards. A chain of resonant plan¬ 
ets, which is a potential outcome of disk-driven migra- 


tion ( Terquem fc Papaloizou|20^ Ogihara fc Ida|2009 ), 
could have similarly interesting consequences on planets 
in the HZ. 

All of these scenarios, which involve coupling be¬ 
tween atmospheric mass loss, tidal evolution, and secular 
planet-planet interactions, probably occur even for plan¬ 
ets outside the HZ. Modeling this coupling c ould be criti¬ 
cal to understanding systems like Kepler-36 (Carter et al. 
20121 , where two planets with semi-major axes differing 
by ^ 10% have a density ratio close to 8, which Lopez fc 


Fortney (2013) show could be the result of starkly differ- 
ent mass loss histories. Future work should investigate 
how mass loss modifies the orbital interactions in multi¬ 
planet systems. 


decreases the IHZ distance. Since evaporated cores are 
more likely at smaller a, a closer-in inner edge to the HZ 
could greatly increase the parameter spa ce a vailable to 
these planets. However, as we argued in §5.2[ HECs are 
li kely to have abundant surface water. 


Luger (fe Barnes (2015) showed that terrestrial planets 
in the HZ of iVl dwarfs may experience long runaway 
greenhouses during their host stars’ extended pre-main 
sequence phases, during which time water loss to space 
can lead to their complete desiccation, rendering them 
uninhabitable. HECs are naturally more robust against 
severe water loss, given that their initially dense H/He 
envelopes can shield the surface and lower atmosphere 
from XUV radiation. Once the envelope is lost, water 
loss from the surface could occur, but by that point the 
stellar XUV flux will be much lower. Future work will 
address the fate of the water on HECs in detail. 


5.9.2. Thermal Evolution 

As mentioned before, we do not consider how the mass 
loss rate affects the planet radius, but instead assume 
the radius instantaneously returns to the “equilibrium” 
value for the given age, mass, and hydrogen fraction. In 
reality, the radius is likely to remai n inflated for som e 
time, particularly for fast mass loss (Lopez et al.||2012). 
Because our radii decrease too quickly, the mass loss 
feedback on the tidal evolution is always negative: mass 
loss always acts to dampen the tidal evolution. How¬ 
ever, accurate modeling of the radius could enable a pos¬ 
itive feedback, in which the effect of the decreasing mass 
(higher da/dt) overpow ers the ef fect of the decreasing 
radius (lower da/dt) in (Bl) and (Cl), leading to faster 
orbital decay. 

Currently, our planet radii are also independent of both 
the insolation and the degree of tidal heating. Under 
the extremely high fluxes and strong tidal forces at early 
times, planets in the HZs of M dwarfs are likely to be 
significantly more inflated than modeled here, resulting 
in faster atmospheric mass loss and likely a higher prob¬ 
ability of complete envelope loss. We also do not model 
radiogenic heating, which could further add to the infla¬ 
tion of the planet. In this sense, our results are conserva¬ 
tive, and HECs may in fact be possible at larger distances 
from their parent stars. A self-consistent thermal evolu¬ 
tion model must be developed to accurately address this 
point. 

We also emphasize that the quasi-static approxima¬ 
tion (see the Appendix) may not be valid for the most 
inflated planets on eccentric orbits above e « 0.4. Given 
the large uncertainty concerning the orbital migration 
due to RLO for planets on eccentric orbits, we urge cau¬ 
tion in interpreting our results quantitatively for highly 
eccentric planets. 


5.9. Other Caveats 
5.9.1. The Habitable Zone 


Recently, Yang et al. (2013) argued that cloud feed¬ 
back on tidally locked planets can greatly increase the 
planetary albedo a nd move the IHZ in by a substantial 
amount. Similarly, 


ets with limited sur: 


Abe et al. (2011) showed that plan- 


ace water are stable against runaway 


greenhouses at greater insolation than Earth, which also 


5.9.3. Tidal Evolution 

Given that the CPL tidal model is accurate only to sec¬ 
ond order in eccentricity, most of our integrations were 
performed under the CTL framework. At low e, however, 
these models predict qualitatively similar results for both 
the orbital migration and the coupling to the mass loss 
processes. At high e, on the other hand, tidal models 
lack observational validation, given the relatively low ec¬ 
centricities of the major bodies in our solar system. We 
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again urge caution in interpreting quantitative results for 
e > 0.5. 

Recall that our values of Qp and Tp are fixed through¬ 
out the integrations and change discontinuously to the 
more dissipative rocky values as soon as the envelope is 
lost; this is certainly an oversimplification. Tho ugh we 
defined Q via its relationship to the phase lag in (19), it 


may also be defined as t he specific dissipation function 
( Goldreich &: Soter||1966 ), defined as 


^ 271 - 


Q-^ = 


2TrEi 


0 Jo 


dE\ , 


(26) 


where Eq is the maximum tidal energy stored in the sys¬ 
tem and the integral is the energy dissipated per cy¬ 
cle. Q is therefore inversely proportional to the dissi¬ 
pation per unit energy stored. Planets with large val¬ 
ues of fn, where the dissipation throughout most of the 
body is small, should therefore have large values of Qp; 
rocky planets with small fn, conversely, have large spe¬ 
cific dissipation rates and correspondingly low values of 
Qp, in agreement with the remarks in (2.4.1 Therefore, 
we would expect that as planets lose mass, Qp should 
decrease, corresponding to an increase in Tp. Modeling 
Qpifn) and Tp{fH), however, is outside the scope of this 
paper, and will be considered in future work. We sim¬ 
ply note that as planets lose mass, the decreasing value 
of Qp/increasing value of Tp should offset the decreasing 
radii, leading to stronger tidal evolution and stronger 
couplings than reported here. 

One must likewise be careful in choosing the planet 
radius that goes into calculating the tidal evolution. This 
radius should be the effective radius of the dissipating 
material, and is thus dependent on the chosen value of Qp 
and Tp; for high Qp/low Tp, the radius should probably be 
that of the envelope. We took this to be Rp, the 20 mbar 
radius shown in the tracks in Figure other choices of 
the tidal radius will lead to different evolutions. 

We only consider planets with zero initial obliquity. 
For nonzero obliquity, the tidal evolution equations must 
be modified and will lead to differences in the evolution; 


see Heller et al. (2011) and Barnes et al. (2013). We also 
ignore the spin-up of the planet as it thermally contracts 
over time, since it should return to the equilibrium spin 
almost instantaneously. The excess angular momentum 
would most likely be absorbed into the orbit, but it is a 
small enough fraction of the orbital angular momentum 
that it can be safely neglected. 

Finally, we note that throughout this paper we set 
tstop = 5 Gyr, based on the age of the solar system. 
For systems older than 5 Gyr, our results should change 
little, since both the tidal and atmospheric evolution will 
have strongly tapered off by this time. Younger systems, 
however, may still be in the throes of tidal decay and 
mass loss processes, and this must be kept in mind when 
searching for HEGs. 


6. CONCLUSIONS 

We have shown that gas-rich mini-Neptunes that mi¬ 
grate early {t < 10 Myr) into the HZs of M dwarfs 
can naturally shed their gaseous envelopes and form gas- 
free Earth-mass planets. Together, Roche lobe overflow 
and hydrodynamic escape can remove up to a few Earth 
masses of hydrogen and helium from these planets, po¬ 


tentially turning them into “habitable evaporated cores” 
(HEGs). This process is most likely for mini-Neptunes 
with solid cores on the order of 1M 0 and up to about 
50% H/He by mass, and can occur around all M dwarfs, 
particularly close to the inner edge fo the HZ. HEGs are 
less likely to form around K and G dwarfs because of 
these stars’ shorter super-luminous pre-main sequence 
phases and shorter XUV saturation timescales. Fur¬ 
thermore, we find that HEGs cannot form from mini- 
Neptunes with core masses greater than about 2 M 0 and 
more than a few percent H/He by mass; thus, massive 
terrestrial super-Earths currently in the HZs of M dwarfs 
have probably always been terrestrial. Ou r results are 


thus similar to those of Lammer et al. 
that planets more massive than 


(2014), who showed 
TTb M 0 typically can¬ 


not lose their accreted nebular gas in the HZs of solar- 
type stars. 

Whether or not a given mini-Neptune forms a HEG 
is highly dependent on the early XUV evolution of the 
host star. In particular, a long XUV saturation timescale 
(Gat ^ 1 Gyr) is needed to fully evaporate the envelopes 
of mini-Neptunes in the HZs of early and mid M dwarfs. 
While a la rge Gat is consist ent with the long activity 
timescales (West et al . 2008) and long spin-down times 
(Pizzolato et al. 2000) of M dwarfs, more observations 
are needed to pin down tsat for these stars. Moreover, 
the relative strength of the X-ray and EUV luminosity 
early on also affects whether or not HEGs can form, 
as this determines whether the atmospheric escape is 
energy-limited or radiation/recombination-limited. In 
the energy-limited regime, which occurs for an X-ray 
dominated flow, the escape is fast and HEGs can form 
throughout most of the HZ of all M dwarfs. In the 
radiation/recombination-limited regime, which applies 
to an EUV dominated flow, HEGs only form close to 
the inner edge of low mass M dwarfs. 

We further find that HEGs can form from planets 
on circular as well as eccentric orbits, though they are 
marginally more likely to have higher e in the long run. 
While there exist feedbacks between atmospheric mass 
loss and tidal evolution, we find that these are only sig¬ 
nificant at e > 0.5; in these cases, whether or not a HEG 
forms can depend just as strongly on the tidal proper¬ 
ties of the planet as on the efficiency of the atmospheric 
escape. 

Many of the Earth-mass terrestrial planets detected 
in the HZs of M dwarfs in the coming years could be 
habitable evaporated cores. These planets should have 
abundant surface water and are likely to be water worlds, 
whose potential for habitability should be investigated 
further. 
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APPENDIX 

MASS LOSS FOR ECCENTRIC ORBITS 
The Quasi-Static Approximation 

In §2.3.6| we argued that for planets on eccentric orbits, the Roche lobe radius may be calculated by replacing a with 
the instantaneous planet-star distance r in ([^, leading to tidally enhanced mass loss close to pericenter. However, this 
is true only if the orbital period is large compared to the dynamical timescale of the pla net, thus allowing suffi cient 
time for the atmosphere to assume the equilibrium shape due to the changing potential. Sepinsky et al. (|2007) call 
this lim it the quasi-static app roximation. Below we show that this approximation is valid for most of our planets. 

As in Sepinsky et al. (2007), we begin by introducing the parameter 




where 


/ = 


(1 - e)3/2 


Porb (1 - e)3/2 

Prot (1 + e)l/2 


(Al) 


(A2) 


is the ratio of the rotational angular velocity to the orbital angular velocity at periastron. Sepinsky et al. (2007) show 
that provided the condition 


Porh 

Tdyn 


> a(e,/) 


(A3) 


holds, the system may be treated as quasi-static. The orbital period is given by Kepler’s law. 


Porh = StT < 


GM* 


(A4) 


while Tdyn, the dynamical timescale of the planet, is 


Tdyn 




(A5) 


The ratio of the two will be smallest for close-in, low-mass planets with large radii. The minimum value in our runs 
occurs for a super-inflated 2M0, 30 R® planet in the IHZ of a 0.08 Mq M dwarf, for which Porb/Tdyn ~ 2.4. 

We must no w co mpar e th is to a(e,/). The equilibrium rotational period for a synchronously-rotating planet is 
obtained from (B6) and (C6) in the CPL and CTL models, respectively, so that we may write a{e) as 




(1 - e)3G ^2,eq 
(l-ge)i/2 n 


(A6) 


This equation is plotted in Figure 12 Note that a only begins to approach Porb/Tdyn for e > 0.6 in the CPL model and 


e > 0.4 in the CTL model, and only for the most inflated planets in the IHZ. We therefore urge caution in interpreting 
results above e > 0.5, where the quasi-static approximation may not hold for some planets. 


The Mass Loss Enhancement Factor 

Since the flux Fxuv is not constant over the course of one orbit, ([^ must be modified slightly: 


M{t) = 


^xuv^xuvAxuv 


4GM„ 


rit? 


1 - 


m) my 


1 -1 


(A7) 


where r{t) is the instantaneous separation between the centers of mass of the star and the planet and where we have 
plugged-in for Fx\jv in terms of Lxuv and made use of (1^. The parameter ^ must also be modified, as the value of 
the Roche radius will also change as the planet’s distance from the star changes during one orbit: 


m = 


Rro 


che 


^xuv 
/ Mp 


V3M* 


1/3 


Rxuv 




(A8) 
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Fig. 12.— The|Sepinsky et al.|||2007|l parameter o as a function of e for a synchronously-rotating planet in the CPL model (black) and 
in the CTL model (gray), i tie dasned. line corresponds to the minimum value of the ratio -Porb/'^dyn across all our runs. Note that for 
6 $ 0-5, the quasi-static approximation is probably valid. 


To avoid confusion, henceforth ^ will denote the original expression for circular orbits (Equations and |^, while ^{t) 
is the time-dependent parameter given by the expression above. The time-average of M over one orbit is 

mt = ^ / M{t)dt (A9) 

27 r/n Jo 

where n is the mean motion. Now, we know that the relationship between n and the eccentric anomaly E is 

nt = E —e sin E (AlO) 

and that 


r{E) = a(l — ecosE), 


so 


dt 1 , , r{E) 

— = - l-ecosE = -. 

dL n an 


Substituting into (A9|, we have 


{M)t = 


2t: I n 


27ra 


M{E)r{E)dE. 


Now, introducing the mass loss rate for e = 0 and Kt\de = 1 from ([^, 

-^xuv^xuv^xuv 


Mo = 


AGMpO? 


it follows from (A71 and (A8 ) that 


M{E) = Moo^ 


r{EY 1 - - 


3 / a 


1 f a 


2? \r{E) J 2^3 \riE) 


-1 


Plugging this into (A13), and using (All), 


M, 




3 1 fa 


{M)t = ^a r(£;)--a+- - r{E) 


27r 

^ [ 

271 Jo 

Mo 

K’ 


2 ^ ' 2 VC 


(1 — ecos A) — ^ + 


-1 


dE 


2 C 2^3(1-e cos £1)2 


-1 


dE 


(All) 

(A12) 


(A13) 


(AM) 


(A15) 


(A16) 
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where we define the eccentric mass loss enhancement factor 

c27r r 




(1 — ecosE) — ^ + 


2C 2^3(1 - ecos£1)2 


dE. 


(A17) 


Note that for e = 0, ATecc reduces to the circular version of Kude, (^6|. Unfortunately, the integral (A16) is not 


analytic. However, if the last term in the integral is small compared to rhe first two, an analytic solution is possible. 
In particular, if we write 


1 


2^3(1 _ ecosEY 


= 1 


(1 — ecosE) -- 

^ ^ 2f 


(A18) 


then we may neglect the last term provided r] 1. Since the term is largest at pericenter {E = 0), we may instead 
require that 


0 < 


1 


2f3(i_e)2(i_ A_e) 


< 1 . 


(A19) 


If this holds, (AI6) simplifies to 


(M). - 


1 — — — e cos E I dE 


-1 


Mn 


^>277 


2^(1 - t) Jo 

Mn : 


I - 


I - ^ 

, 2?, 


cosE 


-1 


dE 


I - ^ 
2{ 




Mn 


1 _ 3 _ _9_ _ „2 
A j 452 e 


1/2- 


Thus, provided condition (AI9) holds (typically for ^ > 10), we may write 




i_3_A_ 2 
C 4^2 ® 


(A20) 


(A21) 


We note that for e = 0 and keeping only first-order terms in the above expression agrees with that of|Erkaev et al.| 

([2007| for C » 1. _ _ 

In Figure 13 we plot the loss rate enhancement factor of Erkaev et al. (2013), l/ATtide (blue lines), along with the 


eccentric version, l/ATecc (exac10 expression in black, approximate expression in red). Recall that there are two distinct 
effects contributing to the extra mass loss for eccentric orbits: the l/Vl — e2 flux enhancement and the smaller Roche 
lobe distance during part of the orbit. In order to better distinguish between these two, we display both l/ATtide 
(dotted blue lines) and I/(Aitide'\/l — e2) (dashed blue lines). Thus, any difference between l/ATtide and I/ATecc in the 
figure is due solely to the Roche lobe effect. 

Let us first compare the curves corresponding to the exact expressions (blue and black). For e = 0, ATecc = Ktide, as 
expected. As the eccentricity increases to 0.25, the flux enhancement effect remains small (i.e., the dashed and dotted 
curves are similar), while the Roche lobe effect begins to become significant (the dashed black curve exceeds the blue 
dashed curve), particularly for low values of A For e ^ 0.5, the effect becomes even more important, leading to an 
enhancement of a factor of several for ^ < 10. 

A particularly important effect is that high e c centri cities will cause the planet to undergo Roche lobe overflow at 
larger values of In the circular jErkaev et al. (2007) model, Roche lobe overflow occurs when ^ = 1 by definition. 
However, it is straightforward to snow that for 


^ — Ccrit — 


1 - e 


(A22) 


the expression in the integral (A17) diverges, resulting in an infinite mass loss rate. This occurs because at pericenter, 
where r = a(l — e), ^(r) = 1. in other words, for ^ < ^crit, the planet will overflow its Roche lobe during at least part 
of its orbit, leading to rapid mass loss. 


3 While we refer to jA17| as the “exact” expression for the (in¬ 
verse of the) mass loss enhancement factor, it is important to re¬ 


member th at it is the eccentric v ersion of the third-order expression 
derived by |Erkaev et al.| (|2007|l and is, in this sense, still an ap- 
proximatiori to tne true ennancement. 
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Fig. 13. — L oss rate enhancemen t factor l/A" as a function of the normalized Roche lobe distance ^ = Rvoche/Rp for 6 = 0, 0.25, 0.5 
(2007i model (.Rtide) is plotted in bl ue, with the flux enhancement factor 1/y/l — (da shed) and without 


Erkaev et al. 


and 0.75. The 
it (dotted). The exact expression 


vprr d erived in this work ( |A16| ) is plotted in black and the analytic approximation ||A20[ ) is plotted in 
red. For low values of condition ( |A19[ | is not satisfied and tne approximation diverges significantly. However, at high eccentricity, as ^ 
increases, the analytic approximation converges to the exact value of Aecc sooner than i^tide- Finally, we note that for sufficiently large 
all curves accounting for the flux enhancement converge to the same value. 


It is important to note that the approximate expression (dashed red curve) converges only for ^ > 10. For smaller 
values of it overestimates the mass loss enhancement significantly and should therefore not be used. However, for 
high eccentricities (see the last panel, for instance), it conve rges to l/itTecc quicker than l/iCtide- Provided condition 
(A19) is satisfied, the approximate analytic expression (A20l does a good job at modeling the actual mass loss. 


Orbits Crossing the Critical Radius 

As we described in §2.3.5[ the escape regime for an EUV-dominated flow may switch from energy-limited to 
radiation/recombination-limited above a certain critical value of the flux, Ecrit, given by 


where 



^ _ ttcxuvAxuv 
GMpi^tide 


(A23) 


(A24) 


and 


B = 7.11 X 10^ g^s” 


Ep 

Ra 


(A25) 


are the coefficients mult iply ing the flux in the energy-limited and radiation/ recombination-limited equations, respec¬ 
tively (Equations and 13). At fixed XUV luminosity, this corresponds to a certain critical orbital radius. 


?"crit = 


Lxuv 

47rE,Ht, 


(A26) 


Planets on low-eccentricity orbits that do not cross rcrit are therefore safely within either the energy-limited (a S> rcHt) 
or radiation/recombination-limited (a rent) regime. In this case, the orbit-averaged mass loss rate is determined 
simply by replacing E in § and ( |T3| ) by its orbit-averaged value, (E), given by (0. 

However, when an orbit is sufficiently eccentric that the atmospheric escape regime switches from energy-limited to 
radiation/recombination-limited over the course of a single orbit, the method outlined above is no longer rigorously 
correct. We must instead integrate the two mass loss rate expressions over the portions of the orbit where they apply. 
Fortunately, as we demonstrate below, these integrals are analytic. This allows us to calculate the time-averaged value 
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of the mass loss rate, {M)t, much as we did in i 2.3.6 


A 


B 


{M)t = pl Fit)dt+- 


lEL 


P 


F 2 {t)dt, 


(A27) 


'RRL 


where P is the orbital period, F{t) is the instantaneous XUV flux, and EL and RRL corr espon d to the energy-limited 
and radiation/recombination-limited portions of the orbit, respectively. By noting as in (A12) that 


an zna 


(A28) 


where E is the eccentric anomaly, we may write 


mt = 


B 

2'Ka 


F^{E)r{E)dE 


r27r—Ec, 


A 

— I E{E)r{E)dE + 

z'Ka ,1 


B 

27ra , 




27r-£;crit 


F^{E)r{E)dE, 


(A29) 


where 


^crit = COS ^ (A30) 

\e ae J 

is the value of the eccentric anomaly when r = Tcrit- The three integrals above follow from the fact that starting 
from pericenter (E = 0), the planet is (by assumption) in the RRL regime up until E = Acrit, switches to EL until 
E = 2tt — Ecrit, and completes the orbit in the RRL regime. By symmetry of the orbit, the first and last integrals are 
identical, so we may simplify: 


(M)t = 


A 

27ra . 


^27r—iJcrit 


F{E)r{E)dE 


— / F^{E)r{E)dE. 
Tro ./n 


Now, noting that E{E) = Lxuv/47rr^(if) and r{E) = a(l — ecosE), we have 

ALxvv 


r2TT-Ecrit 


{M)t = 


dE 


B L 




1 — e cos E 7ra 


XUV 

47r 


dE 


(A31) 


ALxvv dE BEcrit j Axuv 

Acrit 1 —ecosE tto V 47r 


By evaluating the integral, we may finally write 


{M)t = AF 


1-tan 

TT 


-1 


(1 -I- e)tan(Acrit/2) 


BF^ 


(1 - e^)‘^Ecrit 


(A32) 


(A33) 


Note, importantly, that the expressions above are valid only for planets that cross the critical radius during their orbit; 
that is, planets for which the expression a(l — e) < rcut < a(l -I- e) holds. The mass loss rate for planets outside this 
region must be calculated via the method described at the beginning of the section. 

Finally, note that the formalism derived here makes use of A^tide rather than its eccentric version ATecc, which we 
derived in §A.2[ It is in principle possible to account for the tidal enhancement during the energy-limited portion of 
the orbit, but the resulting integral would not be analytic. Moreover, the tidal enhancement due to the eccentricity is 
important primarily near pericenter, where the mass loss is radiation/recombination-limited and therefore independent 
of ATtide- Thus, while the method outlined above may underpredict the mass loss rate in some cases, the effect will be 
small. 
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TIDAL EVOLUTION EXPRESSIONS: CPL 

For zero inclinati on and zero obliquity, the tidal evolution equations for a, e, and uji in a two-body system are 
( Barnes et al.|[2M3 ) 


da 

dt 








147 1 

~20eo,i -I—^^14 ~ 


de 

dt 


ae 

SGM^Mp 


f 2eoy - 



1 

2 


£2,1 + 3e 



duji 

dt 


Z[ 


SM^rl^Rjn 


(dfio.i + [~20eo,i + 49eiy -I- e2,i]) , 


(Bl) 

(B2) 

(B3) 


where the sums are taken over the two bodies. Here, G is the gravitational constant, Ri is the radius of the body 
(planet or star), and Vg^i are the radii of gyration. The parameter is defined 


Z' = 3G'^k2,^Mf{M, + Mg) 


nQi' 


(B4) 


where k 2 ,i is the Love number of degree 2 of the body, which is of order unity and quantifies the contribution of 
the tidal deformation to the total potential, n is the mean motion of the secondary body (i.e., the planet), and Qi are 
the tidal quality factors. The parameters eN,i are the signs of the phase lags (assumed equal in magnitude) of the 
wave on the body, calculated from 


eo,i = sgn(2a;i - 2n) 
ei,i = sgn(2a;i - 3n) 
e 2 ,i = sgn(2a;i - n) 
e 5 ,i = sgn(n) 
es,i = sgn(wj - 2n) 

eg.i = sgn(wj). 


(B5) 


Since shor t-period planets around M dwarfs are likely to be tidally loc ked, one need n ot calculate the planetary spin 
from (B3). Instead, the planet’s rotation rate may be calculated from ( Goldreidi|| 19661: 


=n{l + 9.5e^), 


(B 6 ) 


where n is the mean motion. 


The Typical Case 

Because of the fast rotation of M dwarfs at early times, planets in the HZ are likely to be far outside the corotation 
radius of their parent stars. It is useful at this point to con side r as an example the specific case of a tidally-locked 
planet for which n <C w*. In this case, the stellar phase lags (B5) are all positive and cn,* = +1. For a tidally locked 
planet, e 2 ,p = £ 5 ,^ = £ 9 ,^ = 1, and ei^p = es,p = —1. The parameter eo,p, however, is less straightfo rward to calculate. 
If tida l locking is to be maintained, the average angular acceleration over one period must be zero. |Ferraz-Mello et ^ 


(2008) show that the only self-consistent way to ensure this is if £ 9,2 has a different magnitude than the other lags, 
equal to 


which in our notation corresponds to 


£-o,p — 12e £^ 2 ,p; 


£o,p — +12e . 


n in (B5). 


(B7) 


(B 8 ) 


Note that if e = 0, eo,p = 0, consistent with cOp = 
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In the limit M* ^ Mp and keeping only terms up to order e^, equations (Bl) - (B3) now reduce to 




- 2 iyGM 3 


k2,pRl 

QpMp 


21 

y 




k2,pRp 

QpMp 


a-13/2 


- 13/2 


dt 


3 G^M^k 2 ,,Rl f ^ 2 \ 
32 \ 2 ^ ) 


a 


-9 


(B9) 

(BIO) 

(Bll) 


We need not calc ulate dw^Jdt, since we assume the planet’s spin is instantaneously set to the equilibrium value. The 
first term in (B9| and (BIO) is the orbital effect of the tide raised by the planet on the star. In both equations this 
term is positive, implying that the stellar tide acts to increase both a and e. This makes sense under the assumption 
that the planet is outside of corotation; the stellar bulge therefore leads the planet in the orbit, torquing the planet 
such that it speeds up and migrates outwards. For an eccentric orbit, the strongest impulse occurs at pericenter; since 
the planet must return to that point in the orbit, the pericenter distance is preserved, but the faster orbital speed 
results in a more distant apocenter, and thus higher eccentricity. A similar analysis for the tide raised by the star on 
the planet (the second term in each of the above equations) yields the result that such a tide acts to decrease both a 
and e. 

The evolution of the orbit will therefore depend on the relative magnitudes of the stellar and planetary tides. A 
simple order of magnitude calculation will show that the ratio of the tide generated by the star on the planet to the 
tide generated by the planet on the star is 





(B12) 


For a 10M 0 , 2 R 0 mini-Neptune with Qp = 10"^ orbiting a 0.1 Mq, 0.15R© star with Q* = 10®, this becomes 


^•k—kp 

O^p—kk: 


3 X 10^ 


l + fe2 


(B13) 


which is 1 for all e > 0.01. For an Earth-like planet with Qp = 10^, the ratio is 1 for all e > 0.001. Therefore, 
for planets starting with eccentricities > 0 . 1 , we would expect the planetary tide to dominate the evolution, such that 
the planet’s orbit will shrink. 


TIDAL EVOLUTION EXPRESSIONS: CTL 


The expressions for the evolutions of the orbital parameters are 


( Barnes et al.||2013 | 


da _ 2a? ry ( / 2 (e) /i(e) \ 

dt ~ GM^Mp * V/3^^(e) n /3i5(e )) 

de _ line 7 f /4(e) w* 18 / 3 (e) \ 

dt ~ 2GM^Mp ^ * Vy°(e) n 11 ) 

_ y / / 2 (e) _ / 5 (e) o’A 
dt V/3^^(e) P^{e) n J ’ 


Z, = 3G^k2,^MUM, + 


(Cl) 

(C2) 

(C3) 


(C4) 


where 
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and 


/3(e) 

/i(e) 

/ 2 (e) 

/3(e) 

/4(e) 

/5(e) 


\/l — 


1 + 


31 


1 + 


15 


255 


45 


1 + 


8 

3 2 , 1 4 


185 


16 

, 15 2 15 4 5 

' + +64^ 


l + 3e" + -e^ 

O 



As before, if we assume the planet is tidally locked, its rotation rate is given by 




CTL 

p,eq 


= n 


/2(e) \ 

/3^(e)/5(e)y ■ 


(C5) 


(C 6 ) 


The Typical Case 

As in the CPL mod el, w e exp ect t hat f or th e typical case the tides raised on the planet will dominate the evolution. 
Plugging the result of (C 6 ) into (Cl) and (C2), we find that the second (negative) terms dominate and the orbit should 
therefore shrink and circularize. Proceeding as before, we have 





where 


F 



file) 

/5(e)/i(e) 

/i(e) n 

241 



- 1 


- 1 



for small e. For the same mini-Neptune considered in the 


CPL case (B13), this becomes 




5 X 10® 



(C7) 


(C 8 ) 


(C9) 


which is 3> 1 for all e > O.OOiyw*/ n — 1, implying that the CTL model also predicts a net inward migration due to 
tides. 


RATE OF CHANGE OF THE FLUX 

Conservation of angular momentum requires that for tidal evolution, the rates of change of the eccentricity and the 
semi-major axis must be related through 


1 de 1 da 
e dt 2a dt 


(Dl) 
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The rate of change of the flux due to orbital changes is 
dF _ Lboi d '' ~ ^ -1 

dt 47r dt 

Fhol 


4:TT 

-^bol 1 

47r 


= -F- 
dt 


(a^V^ 

(a^VT^ 


-2 


„ da n -7 , 

2a—vl — 
dt 


2 da 


1 


( 

V2v/1 - e2 

de 


a dt 
4 


(1 


_ g2^3/2 


= -F 


de 

dt 


4 — 5e^ 


e(l - e2)3/2 


For e < 0.7, this simplifies to 


dF ^ -^1 de 
dt e dt' 


which is equivalent to the trivial result 


1 dF 
F dt 


2 da 
a dt 


„ de 
—2e — 
dt 


(D2) 


(D3) 


(D4) 


Thus at lo w e ccentricities, the flux will always increase as the orbit shrinks. However, at very high eccentricities 
(e > 0.8), (D2) predicts that dF/dt is negative when da/dt < 0: the decrease in the flux due to the circularization of 
the orbit overpowers the increase due to the shrinking orbit. 
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